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ABSTRACT 

We investigate the influence of the initial over-densities and masses of proto- 
galaxies on their subsequent evolution (the star formation history in particular), to 
understand whether these key parameters are sufficient to account for the varied prop- 
erties of the galactic populations. By means of fully hydrodynamical N-body simula- 
tions performed with the code EvoL, we produce twelve self-similar models of early- 
type galaxies of different initial masses and over-densities, and follow their evolution 
from the early epochs (detachment from the linear regime and Hubble flow at z > 20) 
down to the stage when mass assembly is complete i.e. z < 1 (in some cases the models 
are calculated up to z — 0. The simulations include radiative cooling, star formation, 
stellar energy feedback, re-ionizing photo-heating background, and chemical enrich- 
ment of the interstellar medium; we do not consider the possible presence of Active 
Nuclei. We find a strong correlation between the initial properties of the proto-haloes 
and their subsequent star formation histories. Massive (M tot ~ 1O 13 M0) haloes ex- 
perience a single, intense burst of star formation (with rates > lO 3 Af /yr) at early 
epochs, consistently with observations, with less pronounced dependence on the initial 
over-density; intermediate mass (Mtot — 10 11 M Q ) haloes have histories that strongly 
depend on their initial over-density, whereas low mass haloes (Mtot — 1O 9 M ) always 
have erratic, bursting like star forming histories, due to the "galactic breathing" phe- 
nomenon. The model galaxies have morphological, structural, and chemical properties 
resembling those of real galaxies, even though some disagreement still occurs, likely a 
consequence of some numerical choices. We conclude that total mass and initial over- 
density drive the star formation histories of early type galaxies. The model galaxies 
belong to the so-called quasi-monolithic (or early hierarchical) scenario in the sense 
the aggregation of lumps of dark and baryonic matter is completed very early on in 
their history. In this picture, nature seems to play the dominant role, whereas nurture 
has a secondary importance. 
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1 INTRODUCTION 

Understanding how early-type galaxies (ETGs) formed and 
evolved along the Hubble time is still a formidable chal- 
lenge. According to the current view of the subject, galax- 
ies probably began to form at z ~ 20 — 50, when Dark 
Matter (DM) haloes containing Baryonic Matter in cos- 
mological proportions gave origin t o the first suffi c iently 
deep gravitational potential wells (|Tegmark et all 1 1991 
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iGao et al.|[2007al lbl). However, it is not yet clear how the 
paradigm A-CDM cosmological model (which prescribes 
a hierarchical, bottom-top formation of structures, imply- 
ing that massive systems assemble their mass later than 
the smaller ones) could be reconciled with the observa- 
tional evidence for large, and re d, galaxies already in 
place at very high redshifts (see e.g. Marchesini et al]|2010l ; 
lHarrison etail l2011al ; iMortlock et alj|201ll h The problem 
is twofold: on one hand, it must be explained how and 
when star formation is quenched in both massive and 
small haloes, which is necessary to reconcile the theoreti- 
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cal prediction with the observed galaxy mass function (see 
e.g. iBundv et all I2006D ; on the other hand, it should be 
also clarified how such massive systems can form at very 
high redshi fts. Semi-analytical models of hierarchical forma- 
tion (e.g., De Lucia et al.ll2006l: De Lucia fc Blaizotl 
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De Lucia et al.ll201ll ) have succeeded in explaining some of 



the observational features of ETGs, but fail in others. In 
particular they overestimate the number of faint objects and 
underpredict the sizes of bright ETGs. In addition to this, 
the hierarchical formation finds it difficult to explain the 
trend shown by the so-called a-enhancement at varying the 
mass of the galaxy. While massive ETGs are enhanced in 
a-element, the low mass ones are not. This requires that 
the history of star formation varies with the galaxy mass in 
a way that can be hardly met by the hierarchical scheme. 
However, since the top ic has been examined in detail by 
l|Tantalo fc Chiosill2002h and no new arguments have been 
made available we will not discuss it in details (but see Sect. 
l3~3l> . 

All the above issues are essentially a reformulation 
of the long-standing dichotomy between the so-called hi- 
erarchical and monolithic mo dels of galaxy form ation, as 
described respective ly e.g. by lEggen et al.l (|l962l ) and by 
IWhite fc Reesl (|l9T8l ^ , possibly including some more sophis- 
tic ated features (e.g . the dry merger scenario, as proposed 
by iBell et ail 12004 ). After the Cold Dark Matter cosmo- 
logical model has got worldwide acceptance, the hierarchi- 
cal model appeared to be the most convincing scenario of 
galaxy formation. In recent times, however, the situation has 
become by far more intrigued than in the past. Evidences 
claiming for a new scenario very close indeed to the classi- 
cal monolithic one have long b een known both observa tion- 
ally and theoretically (see e.g. IChiosi fc Carraroll2002h, but 
recently they hav e become compelling (e.g. iThomasffeoill : 
iGobat et alJhoill). 

iMerlin fc Chiosil (|2006l . l2007n suggested that a series 
of primordial mergers of small stellar subunits (in a sort 
of early hierarchical galaxy formation mechanism), taking 
place at z > 2, could lead to the formation of massive objects 
in the early Universe, resembling the observed systems as 
far as many different aspects (morphology, density profiles, 
metallicity...) are concerned. This picture finds support in 
numerical simulations of the very first stellar populatio ns in 
the early universe, su ch as those bv lBromm et al.1 (|2002h and 
lYoshida et all (|2003l ) . who found that the first virialized ob- 
jects might form at z ~ 30 with a typical mass of 10 6 Mq. In 
this scheme, a first generation of stars in clumps inside pri- 
mordial haloes of DM would first enrich the medium in met- 
als, and also act as the building blocks of larger systems with 
masses up to 10 12 Mq via gravitational clustering. The mode 
of star formation and the efficiency of it is likely driven by 
the mass and mean density of the proto -galaxy as long ago 

^ (120021) and mo re rece ntly 



pointed out bv lChiosi fc Carrarc 
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argued on observational basis by 

A fundamental piece of the puzzle is the role played 
by energy feedback during the formation of the galac- 
tic systems. Although it has long been known (see e.g. 
iMcKee fc Ostrikeijll977h that supernova explosions can in- 
ject large amounts of kinetic and thermal energy into the 
galactic Interstellar Medium (ISM), only more recently the 
importance of this phenomenon as a source of pressure sup- 



port and large scale turbulent motions is being fully recog- 
nized. The effects of the energy injection are still under in- 
vestigation, since two contrasting consequences are expected 
to follow - i.e. the quenching of the star formation activity 
due to the increasing of the local gas pressure (negative feed- 
back), and the indirect enhancement of the stellar activity 
due to the formation of layers of cold and dense gas within 
the turbulent, shocked ISM (positive feedback). However, on 
larger scales bubbles and super-bubbles of hot and thin ma- 
terial are expected to develop due to the action of energy 
injection, expanding in the ISM and giving rise to galactic 
winds. Detailed simulations by ICeverino fc Klypini (|2009T) 
showed that the correct modeling of the energy injection 
from supernova explosions is sufficient to eject a large frac- 
tion of the gaseous content from a galactic disc. How this 
phenomenon is expected to influence the evolution of mas- 
sive spheroids is still to be clarified. On the other hand, it 
is widely accepted that the energy release from Active Nu- 
clei (AGN) should be a fundamental ingredient for quench- 
ing the star formation activity in such systems. However, in 
this study we show how stellar feedback may be sufficient to 
achieve the goal, provided it is modeled in an efficient way. 

Recently, many numerical studies have focused on the 
formation of individual galaxies instead of large sca l e sim - 
ulations. Among the most recent ones, ISales et al.l (|2010h 
produced different gala xy models with diffe rent feedback 
prescriptions; ISommer^ Larsen fc Toft] i|2010l ) simulated a 
cluster and studied i t s evo lution and the galactic popu- 
lations; IStinson et al.l (|2010h used a multi-resolution tech- 
nique to study the form ation of galaxies within high-mass 
haloes; ICroft et al.1 (|2009l ) examined a set of mod els extrapo- 
lated from a large scale cosmological simulation; iNaab et al.l 
(2007) presented three models produced without feedback of 
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any kind; I Cox e t al. (2006) investigated the merger of spi- 
rals; iKobavashil (j200J, 120051 ) calculated a very large num- 
ber of low r esolution chemo -dynamical models of ellipti- 
cal galaxies; iKawatal l| 19991 ) presented models of ellipti- 
cal g alaxies formed out of co smological perturbations; fi- 
nally, [Ch iosi fc Carrard (I2002T ) presented the first chemo- 
dynamical N-body Tree-SPH (NB-TSPH) simulations of 
galaxy formation based on the monolithic scheme, including 
star formation, cooling, energy feedback, chemical enrich- 
ment and galactic winds at varying the initial total mass and 
density of the proto-galaxy, and brought into evidence that a 
continuous transition of the star formation mode from a sin- 
gle initial episode down to a a recurring series of star forming 
events is possible at varying these two parameters. Each of 
these studies addressed to a particular aspect of galaxy for- 
mation. Although the standard hierarchical scheme is fully 
coherent with the A-CDM cosmology, many observational 
hints seem to indicate that the revised monolithic scheme in 
which very early mergers of subunits made of gas and stars 
are taking place yields theoretical models that are able to 
explain a wide range of observational data for ETGs. There- 
fore, it is worth exploring in more details the possibilities of- 
fered by this scheme as a viable and perhaps complementary 
view of the galaxy formation problem. It might be that in 
the A-CDM context while DM haloes aggregate and merge 
building up ever increasing large scale structures, very early 
on baryons collapse in the gravitational potential wells of 
DM mimicking the monolithic like mode. 

Aims of the present study. We intend to explore 
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the possibilities of the monolithic (or early hierarchical) 
mode of galaxy formation we have just outlined and si- 
m ultaneously to brid g e a ga p between the results presented 
by IChiosi fc Carrard |2002l) and those by I Merlin fe Chios! 
(|2006l . |2007| ~ In the first study, the initial conditions for 
the proto-galaxy were not derived from self consistent cos- 
mological simulations, but were created ad-hoc to model 
a spherical lump of DM + BM (in the right proportions 
depending on the adopted cosmological background) sub- 
ject to£ollap_sej_the DM p articles were distributed following 
the|N avarro et al.l fl996) universal profile, with velocities 
derived from the velocity dispersion a(r) for a spherical, 
isotropi c, and collision-less system with the adopted density 
profile l|Binnev fc Tremainel Tl987); while the BM particles 
(in form of gas at the beginning) were distributed homoge- 
neously inside the DM halo with zero velocity field, mimick- 
ing the infall of prim ordial gas into the potential well of DM 
i| White fc Reesll 1978ft . From this stage, star formation in the 
proto-galaxies was allowed to occur. Finally, supposing that 
at each redshift the over-density of proto-galaxies could vary 
within a suitable interval, many models were calculated at 
varying the total mass and initial over-density. 

In the two studies bv lMerlin fc Chiosil (|200fj . 120071 ). the 
initial conditions for DM and BM particles of the proto- 
galaxies were derived from large scale cosmological simula- 
tions, however adapted to our purposes. Since the procedure 
is much similar to the one adopted here, they will be jointly 
described below. H owever, in those papers no attempt was 
made to confirm the lChiosi fc Carrard (120021 ) scenario about 
the mass and/or initial density dependence of the star for- 
mation history of a galaxy. 

With this target in mind, the NB-TSPH code EvoL was 
fully revised, re- written in parallel language and much im- 
proving the input physics. All details about the physical in- 
gredients and the many tests we have p erformed to validate 
the code can be found in lMerlinl (^OOOMMerlin et al.l (|20ld ) 
to whom the reader should refer. The present pape r com - 
pletes the analysis initiated with | Chiosi fc Carrard (I2002T I 
and continued with lMerlin fc Chiosil (|2006l . |2007| V Here, we 
focus on the relation between the initial conditions of the 
host halo and the final properties of the galaxy making use 
of the parallel code EvoL and a consistent cosmological sce- 
nario to set the initial conditions. We investigate the pro- 
cess leading to the formation of isolated galaxies with differ- 
ent total masses and initial over-densities and cast light on 
whether the internal properties of a "real", cosmologically 
consistent galactic halo made of DM and BM are sufficient 
to obtain the variety of structural and physical properties 
which characterize the population of early- type galaxies an d 
to confirm the predictions made bv lChiosi fc Carrard l|2002l) . 
Furthermore, we will examine whether the size predicted by 
the new models better agree with the observational data 
compared to the results obtained by the semi-analytical hi- 
erarchical models we have mentioned above. 

To this aim, we simulate the formation and evolution 
of twelve isolated early-type galaxy models, from the ini- 
tial stages of the non-linear evolution of the proto-galactic 
haloes at very high redshift, down to at least redshift z = 1, 
and beyond in most cases. To avoid any misunderstanding, 
we remind the reader that the evolution of the density per- 
turbations, once exited from the linear regime, toward to the 
stage of proto-galaxy and of this latter to the present, are 



followed by means of NB-TSPH simulations, the only tool 
fully entitled to this task during the non-linear regime. 

The dichotomy between nature (that is the initial prop- 
erties of a proto-galactic system set ab initio by the primor- 
dial cosmological fluctuations) and nurture (that is, the ac- 
tion of environmental factors) in the process of galaxy forma- 
tion is fiercely discussed nowadays. In this context, we claim 
that the initial over-density of an halo should be generally 
considered part of its nature, since it is essentially fixed by 
the cosmological perturbations from which the halo forms, 
rather than part of the nurture processes, as it is sometimes 
viewed. For example, belonging to a group or a cluster of 
galaxies (that virialize at a certain epoch) imply that the 
mean over-density of the region is high, and therefore it is 
likely that galactic haloes belonging to the region had, on 
average, higher initial over-densities than their counterparts 
forming in the field. Of course, isolated but strongly over- 
dense galactic haloes are possible, but they are less probable. 
If the initial cosmological perturbed density field is a linear 
superposition of independent waves, having the same power 
on all scales in scale-free models, then the highest peaks of 
the density field are expected in regions having, on aver- 
age, larger over-densities, rather than being isolated. In this 
picture, the only "real" nurture processes are those relative 
to the action of nearby objects in terms of feedbacks, ram 
pressure, gas stripping, galactic mergers, etc. The argument 
we made here is that these processes, while being of great 
importance in some cases, are not those leading the basic 
mechanisms of galaxy formation and evolution. 

By varying the initial halo over-density of our models 
(in the way we are going to describe below), we try to mimic 
the action of the "natural" over-density. In some way, one 
could consider our "dense" models as systems belonging to 
generally over-dense regions ultimately becoming groups or 
clusters of galaxies, and our "less dense" models as totally 
isolated (field) galaxies. 

All the model galaxies are obtained from the same ini- 
tial proto-halo, however suitably modified in order to sys- 
tematically change its properties from model to model (see 
below). We also run a few more cases with different pre- 
scriptions, to better explore the effects of some key physical 
quantities and assumptions. 

The paper is subdivided as follows. In Sect. [2] we de- 
scribe the numerical methods and the setting of the initial 
conditions. In Sect. Owe analyze the results of the simula- 
tions. In Sect. [4] we discuss the mass density profiles of the 
stars we have obtained for the models. In Sect. [5] we present 
a preliminary discussion of the Mass-Radius relationship for 
ETGs. Finally, in Sect. [6] we draw some final remarks and 
conclusions. 



2 INPUT PHYSICS OF THE NB-TSPH 
MODEL GALAXIES 

2.1 Setting of the initial conditions 

The ideal procedure to derive the initial structure of a proto- 
galaxy would be to start from large scale cosmological NB 
simulations of typically ~500 Mpc on a side in the frame- 
work of a given cosmological model so that the appearance, 
growth and subsequent aggregation of perturbations of all 
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scales can be suitably described, cfr. e .g. th e Millennium 
Simulation described in ISpringel et al.l (120051 ) who studied 
the formation, evolution and clustering of DM haloes in the 
A-CDM cosmology. Of course these simulations require a 
huge number of particles: for instance, the Millennium run 
followed more than 10 billions particles. Considering the 
huge number of DM haloes (proto-galaxies candidates) that 
come into existence, only a small number of particles will 
be used to describe the internal structure of many of them, 
the smaller ones in particular. If this way of proceeding is 
fully satisfactory from the point of view of a large scale cos- 
mological simulation, it is not viable in the case we want 
to study in detail a single object of galactic size. This in- 
deed will have total masses ranging from 10 8 to 1O 13 M0, 
dimensions from say 1 to several tens of kpc, Mdm/Mbm 
in cosmological ratio where DM is by far dominating, and to 
be properly described will require a large number of parti- 
cles (at least 50,000) for each component. Repeated zooming 
in of smaller sub-portions of the whole initial grid is often 
applied when interested to small size objects such as individ- 
ual galaxies compared to the cosmological volume. However, 
this procedure is quite expensive in terms of computational 
resources. Other strategies can be found that have already 
been adopted in some of the studies we referred to. In partic- 
ular, an isolated perturbation can be artificially created and 
evolved without following the evolution of the much larger 
region of space in which it is located. Moreover, since the 
aim of this study is to explore a wide range of masses and 
initial densities, the creation of an ad-hoc set of haloes is 
the most appropriate method to adopt, instead of searching 
in large scale simulations for the particular haloes suited to 
our purposes. To this aim, we proceed as follows. 

i) We assume the A-CDM concordance c osmology, with 
value s inferred from the WMAP-5 data |Hinshaw et al.l 
120091 ): flat geometry, H = 70.1 km/s/Mpc, Q A = 0.721, 
fli, = 0.046 (giving a baryon ratio of ~ 0.1656), as — 0.817, 
and n = 0.96. 

ii) To describe the growth of primordial p erturbations, 
we use COSMICS, the free software written bv lBertschingerl 
i|l995h . This is particularly easy to use and it has been al- 
ready adopted in many studies as the generator of initial 
conditions for NB simulations. However, we are not inter- 
ested in a full cosmological simulation containing perturba- 
tions at all scales, but only in a portion of it containing 
a perturbation with properties (over-density, mass, and di- 
mensions) fixed a priori. The size of this sub-portion is fixed 
in such a way that the wavelength of the perturbation cor- 
responding to the chosen over-density and mass is similar to 
(but suitably smaller than) the size of the sub-portion. Our 
aim is to construct a reference proto-halo containing DM and 
BM particles in the right proportions each of which with its 
own known mass, position and velocity vector. This refer- 
ence proto-halo, instead of being constructed by hand as in 
IChiosi fc Carrarol (|2002l ). is obtained by means of COSMICS 
in the context of an assigned cosmological scenario. Casting 
the problem in a different way, instead of searching within a 
large scale realistic cosmological box the perturbation most 
suited to our purposes, we suppose that a perturbation with 
the desired properties is already there, and derive the po- 
sitions and velocities of all its DM and BM particles from 
a self-consistent, small-size cosmological box tailored to the 
perturbation we have chosen. 



iv) This reference halo will correspond to the largest 
mass of the whole sample of galaxy models we have cal- 
culated. The grid containing our density perturbation has 
the size of I — 9.2 comoving Mpc on a side, and it is pop- 
ulated by a grid of 46 3 particles. The regular positions of 
the particles are perturbed by COSMICS consistently with 
cosmological random gaussian fluctuations. Furthermore, a 
density peak is constrained to form at some very early epoch 
a virialized structure near to the center of the box (COS- 
MICS allows us to specify the properties of the constrained 
density peak to model the density field in the desired way). 
We impose a gaussian spherical over-density with average 
linear density contrast 5p = 3, smoothed over a region of 
radius 3.5 comoving Mpc. We remind the reader that the 
linear density contrast is defined as < p > /pb g — 1 (where 
pb g is the average matter density of the Universe), and it 
can be extrapolated beyond the linear regime; it would be 
equal to 1.86 at the epoch of virialization. Requiring a higher 
value in COSMICS simply bounces back in time the epoch 
of virialization. 

v) Starting from these input assumptions, COSMICS 
returns the initial comoving positions and the initial peculiar 
velocities of the particles at the time in which the highest 
density perturbation is exiting the linear regime (i.e. it has 
p/pbg — 1 = 1). This is an important point to remember. 
COSMICS is used only up to the exit from the linear regime. 
The subsequent evolution of the DM and BM haloes during 
the non-linear regime is followed by NB-TSPH simulations. 

vi) The small cosmological box containing the perturba- 
tion eligible to become a galaxy is our approximation of the 
local Universe, the one subsequently affecting the evolution 
of the perturbation itself. The validity of this approximation 
will be examined in some detail below. To avoid numerical 
problems with sharp edges, we then single out the sphere 
of radius r = 1/2 centered on the center of the box, change 
the particles coordinates from comoving to the proper phys- 
ical values (this is simply achieved by dividing the comov- 
ing value by the initial cosmological expansion parameter 
a — 1/zi — 1, where z; is the initial redshift provided by 
COSMICS at the end of the linear regime), and add a ra- 
dial outward directed velocity component to each particle. 
This component of the velocity is proportional to the ra- 
dial position of the particle and the initial redshift of the 
simulation. This velocity mimics the effects of the outward 
directed Hubble flow, to take the expansion of the Universe 
into account 

vii) To obtain the initial conditions for other model 
proto-haloes (galaxies) with different mass and/or mean ini- 
tial densities, we start from the reference High Density, High 
Mass (HDHM) halo and proceed as follows: 

Varying the mass. If the spatial dimensions of the box are 
changed, COSMICS modifies in a self-consistent fashion the 
size and consequently the mass of the haloes, leaving, how- 

1 Alternatively, the equations of motion could be modified to ac- 
count for the expansion. However, while this is straightforward 
to do in a comoving frame, it is less commonly done in phys- 
ical coordinates. Moreover, virialized structures (like the model 
galaxies after the initial stages of evolution) are no longer sub- 
ject to the expansion of the Universe and should be described by 
static coordinates. E x amples of the m ethod we adopt are given in 
iKatz fc Gurml [|l99ll) , lKawatal j 19991) and in other similar studies. 
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Figure 1. Initial displacement of gas particles in a central slice 
along the XY plane for model IDIM (distances are in proper 
Mpc). The color code indicates the mass density [g/cm 3 ]. All the 
other models have similar initial conditions. 

ever, the mean density unchanged. Playing with this, we 
generate a number of proto-haloes with the same mean den- 
sity but smaller mass and consequently radius. Considering 
the extremely wide range of galaxy masses observed in the 
Universe, we generated proto-haloes of different mass whose 
value differ by a factor of ~ 64 (and the corresponding radii 
by factors of 4) passing from one to another. The initial red- 
shifts of models with equal initial density but different mass 
are different from one another, because the exit time from 
the linear regime is a function of the mass of the perturba- 
tion. 

Varying the density. The procedure to generate haloes with 
the same mass but different initial mean density is more 
complicate. Before perturbing the positions (and velocities) 
of the particles in the grid, we artificially decrease the den- 
sities at each grid point of the grid, thus varying the average 
density but conserving the total mass. Since the initial con- 
ditions are produced at the moment in which the highest 
perturbation peak exits the linear regime, the initial spatial 
configuration is still a superposition of independent plane 
waves with different wave-numbers and random phases, i.e. 
a gaussian field with S(q) = ^2^S^exp(ik ■ q), where k is 
the wave-number and q is the Lagrangian coordinate corre- 
sponding to the unperturbed comoving position of a mass 
element. Dividing the local (over-)densities 8(q) by a con- 
stant factor / returns a self-similar gaussian field in which 
the Fourier coefficients 5t is divided by the same factor, im- 
plying that the variance a 2 oc Si of the perturbation is 
reduced by a factor f 2 . COSMICS subsequently computes 
the effective displacement and velocities of the particles via 
the Zel'dovich approximation, 

x{q) = q + D+d(q), 
v(q) = D+d(q), 

(1) 

where the displacement field d(q) is given by the relation 
V • d = Sp ^ , and D+ is the cosmic growth factor at 

the initial redshift, which is a function of the cosmological 
model. We avoid the simplest choice of imposing a less pro- 
nounced density peak (i.e., a lower value of the initial density 



contrast 8 in the constrained field), because this would re- 
duce the density of the constraint without modifying the ran- 
dom perturbations, thus resulting in a very different global 
field with more spectral power on the small scales. To ob- 
tain intermediate density haloes, the mean density of the 
reference halo is decreased by a factor of 15%. The mean 
density of these latter is then decreased by another 15% to 
obtain low density haloes. Therefore, these latter have the 
an initial mean density which is 72.25% of the value of the 
initial reference halo. Finally, a set of haloes with very low 
initial over-density was obtained reducing the mean density 
by 50%. As already pointed out, the different densities can 
be considered as a way to mimic the possible environments 
in which a proto-halo begins to form. 

viii) At this stage, a minimal amount of solid-body ro- 
tation is also added to all the particles in the proto-haloes. 
The solid-body rotation is represented by the spin parameter 

A- J|£|1/2 (2) 

where J is the angular momentum, E is the initial binding 
energy, and M the total ma ss of the syste m. Typical values of 
A ranee from 0.02 and 0.08 (|Whitell 19841) . which corresponds 
to angular velocities of the order of fractions of a complete 
rotatio n over time-scales as long as about ten free-fall time- 
scales (|Carraro et al.lll998l ). We adopt A = 0.02. The choice 
is motivated by the fact that our model galaxies are meant to 
represent ETGs. This issue will be also shortly commented 
in general remarks of this section. 

ix) Finally, each particle in each proto-halo is split into 
a DM and a BM particle (gas), this latter with a small dis- 
placement from its original position to avoid numerical di- 
vergences. Indicating with mo the mass of the original par- 
ticle and with fsD — Pbm/pdm — 0.1656 the cosmologi- 
cal ratio between the BM and DM densities, we assign to 
each gas particle the mass m gas — fso x mo- Consequently, 
moM = (1 — 0.1656) x mo. The final HDHM proto-halo con- 
sists of ~ 58,000 DM particles plus an equal number of gas 
particles. 

Coding the models. With this method, we have obtained 
twelve initial haloes. As the reference HDHM model, each 
one is identified by a string of four or five letters, of which 
the first two (or three) refer to its density and the last two 
to its mass(e.g., IDLM is for intermediate density, low mass; 
VLDHM stands for very low density, high mass, etc.). The 
key initial value of a few important parameters for all the 
twelve haloes are listed in Table [1] whereas a sketch of the 
density contrast inside the initial sphere projected on the 
XY plane for the case of the model IDIM is shown in Fig. 
[1] where a number of lumps of matter are clearly evident. 
The number of particles in all models is similar to that of 
the HDHM case, with some minor differences due to the dif- 
ferent displacement of the particles in the outskirts of the 
retailed sphere. 

General remarks. We point out that all the proto-haloes 
we have chosen are fully consistent with the underlying cos- 
mological background. The e xpected DM halo mass fu nc- 
tion follows the la ws by, e.g., IPress fc Schechterl (| 19741) or 
ISheth et al.l l|200ll ). and/or the analyti cal fits of numeri- 
cal simulations (e.g. IWarren et ai]|2006l ). In this picture, a 
M ~ 10 13 M Q halo is becoming "typical" on a ~ 10 Mpc 
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scale (the typical linear size of a large galaxy cluster) at 
redshifts below z ~ 5, which is compatible with the collaps- 
ing redshift of our HDHM model (see below). Such massive 
haloes are also expected not to be rare on ten times larger 
scales at z ~ 10. This can be easily checked considering the 
halo growth functions of DM haloes for the concordance A- 
CDM cosmology: the expected number of haloes with mass 
~ 1O 13 M within a volume of 1 (Mpc /fe) 3 at z ~ 10 is 
n ~ 10~ 8 (see Fig.2 of lLukic et ai1l2007l '). 

The main criticism that could be made to the above 
procedure is that the size of the simulation backing the per- 
turbation in question is too small (9.3 comoving Mpc on a 
side for the HDHM box, and consequently smaller in the 
other cases). It is long known that the simulation size de- 
termines the maximum perturbation wavelength. If the long 
wavelengths are dropped out, the strength of subsequent 
clustering is reduced, but at the same time the number den- 
sity of intermediate mass haloes (with tota l mass of the or- 
der o f 10 13 ft -1 M Q in our case) is enhanced l|Power fc Knebj 
2006). While this may be a problem in real cosmological 
simulations, in our case it is not so, because we are in- 
terested in objects with the mass a nd size of galaxies and 
not of galaxy clusters. Furthermore. IPower fc Knebei (|2006t ) 
show that truncation of the initial power spectrum (i.e., the 
small size) of the simulation has little impact on the inter- 
nal properties of the haloes. However, truncation lead to 
spin parameters that are 15% lower than usual. We have in 
a sense avoided the whole spin problem just by taking the 
lowest value for the spin parameter (A = 0.02) suggested by 
l|Whitelll984l ). a reasonable approximation for slowly rotat- 
ing systems such as ETGs. A final remark could be made 
related to the fact that adopting the sphere of radius r = 1/2 
(about the size of the Local Group) as the starting proto- 
halo we implicitly neglect the possibility that late infall of 
nearby haloes that initially were outside the sphere may oc- 
cur. This would inhibit late refuelling of gas (and stars and 
DM) to our system that is equivalent to say that we inhibit 
late mergers likely accompanied by star formation. This in- 
deed is not the aim of the present study, which intends to 
explore the modality of star formation in alternative to the 
very popular hierarchical scheme. To conclude, the scheme 
we propose to derive the initial conditions for our model 
galaxies is not in conflict with the cosmological paradigm 
and it is easy to implement in NB-TSPH simulations. 



l20ld for the definition), an artificial thermal conduction to 
smooth out contact surface discontinuities, and finally the 
variable adaptive softening lengths with consistent correct- 
ing terms in the equations of motion. In the present models 
we make the following assumptions and/or simplifications: 

(i) We switch off the artificial thermal conduction, be- 
cause its effects in combination with radiative cooling are 
not yet fully understood and tested, and may lead to un- 
physical behaviors, with particles exchanging heat because 
of numerical conduction and simultaneously losing energy 
due to physical cooling. The issue deserves an accurate anal- 
ysis to be planned for a forthcoming study. 

(ii) We adopt variable softening lengths; this is an im- 
portant difference with respect to other similar calculations, 
in which a consta nt softening le ngth has been adopted. The 
only exception is iNelsonl (|2006h even if in that case no cor- 
recting terms in the equations of motion were considered. In 
the present models, we let the softening, e, and smoothing, 
h, lengths vary freely. In a few companion models calculated 
for the sake of comparison, we adopt for the two parameters 
the minimum threshold values e m in = 100 and h m in = 100 
pc. Such low limits would provide sufficient resolution in all 
the cases of interest. However, these side models show that 
the results are very sensitive to these parameters, even if all 
the models are essentially similar in their global aspects (see 
Sect. 13.51 for a discussion on this issue). 

(iii) We include a density-dependent pressure limit to 
avoid artificial clumping of poorly re solved gas. The method 
is sim ilar to the one presented by iRobertson fc Kravtsovl 
(2008). In practice, at each time-step all SPH particles are 
used to compute the local Jeans mass, which is a function 
of density and temperature, i.e. 

_ 47T / / 7T \ 3 

Wljeans — -jj- Pgas [ Cs y ~Q~p J ' ' ' 

where c s — ^7(7 — l)u is the local sound speed, p gas is 
the density of the SPH particle, and ptot is the total local 
density (i.e., gas plus DM and/or stars jfj. Then, the internal 
energy of each particle (limited to the computation of the 
hydrodynamical forces) is chosen as the maximum between 
the real value and the "effective" value, which is defined as 
the energy that the gas sphere should have to be stabilized 
against collapse: 



2.2 The NB-TSPH code EvoL 

EvoL is a Lagrangian N-Body parallel code designed to 
study the evolution of astrophysical systems on any spa- 
tial scale, from large cosmological volumes to sub-galactic 
regions. It is based on the standard Tree Al gorithm for the 
descri ption of the gravitational interactions iBarnes fc Hutl 
(1986), and on a modern version of the Smoothed Particle 
Hydrodynamics method to model gas dynamics. Its basic 
features concerning the treatment of gravitational i nterac- 
tion a nd h ydrodynamics are d escribed in details in iMerlinl 
l|2009h and lMerlin et all (|2010T ); here, we briefly recall some 
of its characteristics. 

EvoL includes some peculiar features, such as the la- 
grangian formulation of the SPH method that includes the 
so-called V/i terms in the equations of motion (taking into 
account variation of the smoothing length, see lMerlin et alj 



u eff = N jeans x — —p tot 5 , (4) 

7(7-1) \AlV2p gas J 

where m is the particle mass, Nsph is the number of par- 
ticles required to resolve a region (indicatively, this number 
can be considered equal to the typical number of neighbours, 
i.e. about 60 in our models, or a multiple of this value), 7 is 
the adiabatic index, and Nj ea ns i s a free parameter > 1. In 
these models we use Nj eans = 15 (Merlin 2009; Merli n et al.l 
l2012l .in preparation) . Note that u e f / is only used to compute 
the mechanical acceleration due to the internal gas pressure. 

2 This expression for the Jeans mass of a gas sphere in presence 
of a non collisional component can be obtained from a pertur- 
bative analysis of the stab ility of a pressurized cloud, see e.g. 
ITittlev fc CouchmanI l|2000h for reference. 
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Table 1. Initial parameters for the twelve models of proto-galaxies. Left to right: total masses, corresponding (initial) gaseous masses, 
mean halo over-densities (<5p — 1) at redshift z = 30, initial redshift, halo initial proper physical radius, mass of a gas particle. 



Model 


Mtot [M ] 




< 5p — 1 > 2 =30 




Tini [kpc] 


rrigas [Af©] 


HDHM 


1.75 x 10" 


2.90 x 10 12 


0.39 


46.34 


97.17 


4.97 X 10 Y 


IDHM 


1.75 x 10 13 


2.90 x 10 12 


0.30 


39.24 


114.31 


4.97 x 10 7 


LDHM 


1.75 x lO 13 


2.90 x 10 12 


0.23 


33.20 


134.49 


4.97 x 10 y 


VLDHM 


1.75 x 10 13 


2.90 x 10 12 


/ 


22.67 


194.34 


4.97 x W 


HDMM 


2.69 x 10 11 


4.45 x 10 IU 


0.46 


53.79 


20.99 


7.79 x 10 5 


IDIM 


2.69 x 10 11 


4.45 x 10 IU 


0.33 


45.57 


24.69 


7.79 x 10 5 


LDIM 


2.69 x 10 11 


4.45 x 10 10 


0.25 


38.59 


29.05 


7.79 x 10 5 


VLDIM 


2.69 x 10 11 


4.45 x 10 IU 


/ 


26.37 


41.98 


7.79 x 10 5 


HDLM 


4.17 x 10 y 


6.91 x 10 s 


0.54 


63.23 


4.48 


1.22 x 10 4 


IDLM 


4.17 x 10 9 


6.91 x 10 s 


0.39 


53.60 


5.27 


1.22 x 10 4 


LDLM 


4.17 x 10 y 


6.91 x 10 s 


0.29 


45.40 


6.20 


1.22 x 10 4 


VLDLM 


4.17 x 10 9 


6.91 x 10 8 


0.16 


31.11 


8.96 


1.22 x 10 4 



In the expressions for the variation of internal energy, dissi- 
pation, and cooling, the real value of u is adopted. In this 
way, the temperature of gas particles can reach low values, 
without reaching high densities, thus creating regions of cold 
material without spurious numerical clumping of particles. 

(iv) We run the simulations in physical coordinates 
since at present we are not interested in large scale mo- 
tions. Moreover, this choice rules out the need for periodic 
boundary conditions, which would require a large amount of 
additional CPU-time. 

(v) Our initial system is confined within the initial 
sphere oil/2 radius and it is surrounded by empty space. In 
other words the evolution of our model galaxy occurs in void. 
This simplification of the problem has some non-negligible 
drawbacks. First of all, the infall of material towards to in- 
nermost regions of the galaxy sooner or later stops because 
material outside the sphere cannot flow in. This rules out 
the possibility of late mergers and/or late infall of signifi- 
cant amounts of gas. Important point to keep in mind, but 
which does not invalidate the results of this study as long 
as late merger, matter acquisition are not of interest. Sec- 
ondly, the feedback-heated gas leaving the galaxy will escape 
freely without experiencing drag forces. Thirdly, the density 
of the system cannot be consistently computed and the mod- 
els evolve as small closed Universes rather than over-dense 
regions within a flat Universe, so the redshifts are only in- 
dicative of the real evolutionary time-line. 

(vi) Finally, we include radiative cooling, star forma- 
ti on, and stella r feedback. As they are described in detail 
in iMerlinl (|2009l ). in the sections below we limit ourselves to 
summarize their key aspects. 



2.3 Radiative cooling 

The cooling functi ons for atom i c rad iative processes are 
those described in lCarraro et al.l (119981) . which i n turn are 
based on those elaborated bv lchiosi et ail (|l998T ). In brief, 
for temperatures greater t han 10 4 K they lean on the 
ISutherland fc Dopital {l993) tabulations for a plasma un- 
der equilibrium conditions and metal abundances \og\Z / Zq\ 
= -10 (no metals), -3, -2, -1.5, -1, -0.5, (solar) and 
0.5. The chemical properties of gas particles are followed 
in details, tracking the abundances of Fe, O, Mg, Si and 
C, and the global metallicity Z. For temperatures in the 
range 100 <T< 10 4 the dominant source of cooling is the 



H2 molecule becoming rotationally and/or vibrationally ex- 
cited through a collision with an H atom or another H2 
molecule and decaying through radiative emission; the data 
in use have been derived from the analytical exp r ession s 
of iHollenbach fc McKeel lll979l) and iTeemark et all | | 19971 ). 
amalgamated together by Chiosi et alj (1 19981 ). We point out 
that in the present simulations it is assumed that all hy- 
drogen becomes molecular below 10 4 K, so the cooling rate 
in this range of temperatures may be somewhat overesti- 
mated. However, cooling from collisions of neutral hydrogen 
atoms with electrons from ionized metals (when present) is 
neglected. This may somewhat compensate for the excess 
of cooling due to molecular hydrogen. In any case this is a 
(marginally) weak point of our current treatment of cooling 

that should be improved. 

Finall y, for temperatures lo wer than 100 K. IChiosi et alj 

ll 19981) and Carraro et al.l (Il998l). starting from the studies of 
ICaimmi fc Seccd (ll98riTandlTheis et all (Il992l ), incorporate 
the results of IHollenbach fc McKeel i| 19791 ) and [Hollenbachl 
( 1988) for CO molecule as the dominant coolant. The follow- 
ing analytical relation in which the mean fractionary abun- 
dance of CO is given as a function of [Fe/H] , is found to fairly 
represent the normalized cooling rate (i.e. Kco/n 2 with n 
the number density of particles) 



^=1.6X 1(r 2 9 10 ([^/H]-1.699)T° 5 (5) 



in erg/s/cm~ 3 . For all other details see lCarraro et al1l|l998h . 

We also include the contribution by inverse Compton 
cooling 



Aeon 



5.41 x 10~ 36 x e [T - T c mb {1 + z)](l + z) 



(6) 

in erg/s/cm 3 , with the usual meaning of symbols (see e.g. 
Ilkeuchi fc Ostrikerlll986l ). 

The radiative cooling is computed separately from the 
SPH equation of energy conservation, due to its very short 
time-scales. During a dynamical time-step, the density and 
the metallicity of gas particles are kept fixed, as well as the 
mechanical heating rate. On the contrary, the temperature 
is let vary under the action of cooling, and the new cooling 
rates are simultaneously obtained as a function of the new 
temperature. The process is iterated until the end of the 
time-step. 
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2.4 Star formation 

Star formation is modeled by m eans of a stocha s tic me thod, 
similar to that introduced by IChurches etHI (|200lh and 
iLia et all (|2002h . First, only gas particles belonging to con- 
vergent flows (i.e. V ■ v < 0) and denser than a suitable 
threshold p» are considered eligible to form stars. In our 
models we adopt p* = 5 x 10~ 25 g/cm 3 . No restriction on 
the temperature is imposed; this is motivated by the fact 
that thermal instabilities can produce star forming sites even 
within high temperature gas. 

If a gas particle satisfies these criteria, it is assumed 
to form stars at the rate dp*/dt = e-SFPg/tff, where tff ~ 
0.5/^/Gptot is the free-fall time, ptot is the local total mass 
density (DM plus BM), and esF the dimensionless efficiency 
of the star formation process. This means that a gas particle 
is expected to transform a fraction esF of its mass into stars 
over its free-fall time scale. However, a stochastic description 
of the star forming process is adopted to avoid the creation 
of exceedingly large numbers of star particles. Thus, gas par- 
ticles undergo a Monte Carlo selection to check whether or 
not they will actually form stars, in such a case they are 
inst antaneously turned into collis ion-less star particles (see 
also ILia et al.l 12002 ; iMer lin 2009). To this aim, a random 
number r £ [0, 1] is drawn and compared to the probability 
that the particle is actually forming stars. An obvious choice 
would be to use the probability P — esF/tff x At, where 
At is the dynamical time-step. However, the evaluation of 
the real probability must take into account that the dynam- 
ical time-steps are generally much shorter than the free-fall 
times of particles. This implies that within a free-fall time 
scale a number of star forming events and corresponding 
random draws are possible depending on At, and therefore 
so does the global probability for a single gas particle to be 
turned into stars. If the global probability over tff is P, the 
stochastic process of random selection must be corrected by 
considering that the probability summed over the number 
of draws n ~ At /tff must be equal to the probability of a 
single draw within the whole free-fall time. First we take as 
the probability within a single time-step the value p ~ P/n, 
then we consider as the probability of a successful event over 
n draws the quantity p x J^" =1 (l—p) !_1 < P. The correcting 
multiplicative term is 



1 p x£"(l -£)<-i- yj 

Therefore, the final probability during the time-step is / x p. 
If r < fp, then the gas particle is instantaneously turned into 
a stellar particle. 

In the present models, we assume e$F = 1- Empirical 
estimates of the efficiency of star formation based on obser- 
vational data of star forming events insi de molecular clouds 
in the local vicinity yie ld e,sF — 0.025 l|Lada fc Laddl2003l ; 
iKrumholz fc Tanll2007l ). In this case the gas density and 
free-fall time scale of are those of the cold molecular clouds. 
However, in the case of the large scale star formation mech- 
anism in a galaxy, the above estimate may not correspond 
to reality. As a matter of facts, considering that a typical 
galaxy with 10 11 M@ mass in stars has to disposal a time 
scale of about 13 Gyr to build up its stellar content, the 
mean estimate of the star formation efficiency is closer to 



0.1 rather than 0.02. The new estimate would increase by 
nearly a factor of ten if the time to disposal to form stars 
is much shorter than 13 Gyr, say 1 to 2 Gyr. In addition 
to this, it is not known whether esF should be maintained 
constant during different cosmic epochs due to the differ- 
ent mechanisms forming stars at high redshifts and in metal 
poor regions with respect to local molecular clouds. Finally, 
it is not known whether the efficiency of star formation is 
the same in all galaxies independently of their mass (either 
total or baryonic). For all these reasons, we prefer to adopt 
here esF = 1, which means that the process of star for- 
mation occurs at 100% efficiency (meaning that all gaseous 
particles satisfying the SF criteria will turn their total mass 
into stars within a free-fall time scale). Another consider- 
ation has to be made here. Several numerical simulations 
calculated with different values of esF clarify that within a 
certain range of values, the star formation histories of the 
models with small values of esF are apparently quite sim- 
ilar to those with esF = 1, the only difference being that 
with small values of esF the simulations become much more 
time-consuming. The reason is that dense and cold clumps 
of matter, interacting with hotter material heated up by 
close SN explosions, require extremely small time-steps. In 
contrast, if the gas can easily form stars, this critical situ- 
ation is avoided, and the simulations proceed much faster. 
The weak impact of esF on the final SF history can be at- 
tributed to the self-regulation cycle between star formation 
and energy feedback. A high SF efficiency implies a strong 
and sudden burst of stellar activity; but young stars soon 
pressurize their surroundings via energy feedback (see be- 
low) , quenching the formation of new stars. However, if the 
gas is sufficiently dense, radiative cooling is effective and 
further star formation can soon take place (positive feed- 
back). With a low efficiency, less stars form. Their heating 
is consequently lower, and more stars can soon form before 
the feedback halts their formation. This ultimately leads to 
the same situation as in the previous case, with perhaps 
the minor consequence of a delayed enrichment in heavy 
elements of the medium. Clearly, other parameters play a 
more important role; for example, the density threshold p, 
and the efficiency of feedback. Therefore the choice esF — 1 
and fast calculations likely, is the best compromise. How- 
ever, one should always keep in mind that adopting esF — 1 
may have other consequences on the dynamical evolution of 
the systems, favouring the collision- less collapse instead of 
the dissipative one (see Sect. [6}. 

2.5 Energy feedback 

When a gas particle is turned into a star particle, it can be 
considered to represent a Single Stellar Population (SSP) 
made of many real stars. It starts re-fueling the ISM with 
heavy chemical elements and energy, mainly because of 
winds from young massive stars and Supernova (SN) ex- 
plosions. 

Winds. We consider two regimes for the energy injec- 
tion by winds. When the SSP is young, the main source 
of energy are the young massive stars. Their winds have 
very high speeds, typically from 1000 to 3000 km/s. Here 
we adopt a constant velocity of vyso = 1500 km/s. The 
kinetic energy of the winds is assumed to be thermalized 
and released within the surrounding medium with an effi- 
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ciency of 20% (see e.g. iDvson fe Williamdll997h . In the late 
stages of the SSP evolution, slow velocity winds from low 
mass stars become predominant; in these case we assume 
the wind speed of voso = 60 km/s (typical of RGB stars), 
and the same thermalization efficiency. 

SN explosions. Each real SN explosion is expected 
to deposit some 10 51 kinetic ergs in a very small region, 
and on a short time scale. However, most of this energy is 
soon radiated away, and only a small fracti on of it is subse- 
quent ly thermalized. We use the results bv lThornton et all 
to obtain an analytic approximation of the fraction 
of the initial energy which becomes available at the end of 
the expanding phase of the SN bubble. However, in the real 
Universe SN explosions take place in already shock-heated 
regions, with temperatures raised up to some 10 6 K, because 
of the action of a photo-ionizing flux from massive young 
stars, stellar winds, and previous SN explosions. All this is 
not taken into consideration in Thornton's study, so a large 
uncertainty still affects the description of the whole process . 
However there are a few exceptions like lCho fc Kang (|2008l ) 
that should be taken into account. A single star particle 
represents an entire SSP, so a large number of SN explo- 
sions are expected to take place within a single star par- 
ticle over a rather long time scale. The number of SN ex- 
plosions (attention must be paid to release the energy in 
discrete bursts consistently with the number of real explo- 
sions) , as well as the amount of gas released by a SSP of given 
age and mass during a time-step and its chemical composi- 
tion a re computed with the technique described in lLia et alj 
(2002), using th eoretical SSPs (in EvoL, the Padova tracks 
are adopted, see iGreggio fc Renzinilll98a) , and adopting the 
initial m ass funct i on fo r the real stars in the SSP (star par- 
ticle) of iKroupal l|200ll ). We refer the reader to the above 
papers for all the details. It is w orth clarifying t hat we do 
not adopt stochastic approach of iLia et all (|2002h to model 
SN explosions. At each time-step they decided the fate of 
each star particle by means of the Monte-Carlo method; 
randomly chosen star particles were instantaneously turned 
back into gas, while the remaining ones continued their life 
as star particles. On the contrary, in our description each 
SSP continuously releases gas and energy at each time step, 
becoming a hybrid particle (in the sense described below). 

Finally, the total budget of thermal energy produced 
by stellar winds and SN explosions is given continuously 
to a sort of gas reservoir assigned to star particles, ulti- 
mately formed by the gaseous material ejected from the SSP. 
This reservoir, due to its high temperature, acts a piston on 
the neighboring particles, effectively injecting kinetic energy 
into the surrounding interstellar medium. In the next Sec- 
tion we clarify some aspects of this method. 

2.6 Fate of the gas ejected by SSPs 

The gas ejected by SSPs, heated up by the energy budget 
from winds and SN explosions and enriched in chemical el- 
ements, interacts with the surrounding medium. This is a 
crucial issue in numerical simulations of galaxy formation. 
Since it is not feasible to continuously form new gas parti- 
cles from the stellar particles, it is necessary to develop a 
consistent, yet simple method to model this process. 

There are two main sources of uncertainty. First, one 
has to decide how to deal with the hot gas soon after its 




Figure 2. Snapshot showing the displacement of gas particles on 
a central slice on the XY plane (model IDIM, z = 7.6; distances 
are in proper Mpc). Features caused by energy injections from 
stars (the clumped brown particles near the origin of coordinates) 
such as expanding bubbles and shock fronts as well as collapsing 
cosmological filaments are clearly visible. 

ejection from the SSPs (this gas is essentially composed by 
matter ejected by stars in form of winds and/or SN explo- 
sions). How should it be let cool? How to compute its den- 
sity? It has long been known that distributing the feedback 
energy among neighboring particles and then letting them 
normally evolve among the other cold ones yields poor re- 
sults and very low feedback efficiencies (Katz 1992) : the en- 
ergy is soon radiated away and, in critical situations, hot 
particles escape from the region due to excessive hydrody- 
namical friction with their neighbors that is not physically 
grounded. Also, distributing the gas released by the SSPs 
to nearby SPH particles at each time-step has proven to 
be inefficient, since the large amounts of energy released lo- 
cally are smoothed out too strongly. To cope with all this, 
we model we proceed as follows. We suppose that initially 
the ejected gas remains stuck to the parent star particles. 
This assumption sounds reasonable in very early stages of 
gas ejection, but surely becomes unphysical as time goes by. 
Sooner or later, the gas ejected by SSPs should move freely 
and mix with its nearby particles. Therefore, the issue is 
how and when the gas should be let evolve "normally" , as a 
standard SPH gas particle - i.e., when it should dynamically 
decouple from its parent star particle. 

To this aim, we split the evolution of the ejected gas 
into two temporal steps. The gas is kept locked to the par- 
ent SSP during the first stages of its evolution, namely from 
its first appearance in the SSP to the age of t e j ec t — 8 x I0 7 
years. This is the typical lifetime of a 6 Mg star with so- 
lar co mposition and convective overshooting (|Bertelli et al.l 
1 1994) . In this type of stellar models, the 6 Mq star is also 
the lower limit for the occurrence of Type II Supernovae. 
The value of t e j ect could be lowered by a factor of 2 (and 
the mass limit accordingly increased to about 8 — 9 Mq ) to 
be more conservative, however without changing the overall 
picture we are proposing. During this first stage, the ejected 
gas does not contribute to the acceleration of the parent 
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Figure 4. Star formation histories of two low mass models with 
(left) and without (right) energy feedback from the stars. 



Figure 3. Phase diagram p (g/cm 3 ) vs. T (K) for the IDLM 
model at z = 8.5, with (left) and without (right) the inclusion 
of a photo-ionizing background radiation at z = 9. The lines are 
the iso-levels of the logarithm of the number of particles at each 
point of the phase diagram, as indicated by the color scale. See 
text for details. 

particle, but acts as a piston exerting pressure on the neigh- 
bouring particles because of its high internal energy. Its den- 
sity is computed normally with the SPH method, and ra- 
diative cooling is let occur normally during this phase. In 
this way, this "almost-zero-mass" SPH particle, locked to a 
collision-less star particle, can effectively act as a source of 
feedback, pressurizing the region and giving rise to complex 
structures resembling observed bubbles and super-bubbles 
as shown in Fig. [2] The me chanism is similar to the one 
proposed bv lPelupessvl ([2005), who used zero-mass pressure 
particles to model the stellar feedback. However, instead of 
introducing zero-mass particles, in EvoL the real gaseous 
mass ejected from SSPs during the current time-step (al- 
ways much smaller than the typical mass of a SPH parti- 
cle) is computed and added to a gaseous reservoir locked 
to the star particles; these "hybrid" particles interact hy- 
drodynamically with the surrounding gas particles via the 
pressure exerted by their gaseous component. 

When tej ect is elapsed, the hybrid star/gas particle is let 
free to collect a sufficient amount of gas (ejected by nearby 
SSPs in their first stage too) from other neighboring hybrid 
particles to form a whole new gas particle. When the mass 
of the collected gas is equal or comparable to the mass of 
a standard SPH particle (in the present models we allow 
~ 50% tolerance), then a new gaseous particle is spawned, 
with position, velocity and temperature obtained from the 
weighed mean over all parent star particles that contributed 
to create it with their gas supply (the mass of gas in these 
hybrid particles is accordingly decreased by the amount of 
gas they have given away). From now on, the new particle is 
let free to evolve normally. This method has proven to give 
satisfactory results, at the expense of a slight increase of the 
total number of particles during the run. 

2.7 Chemical enrichment 

Stellar gaseous ejecta contain heavy elements, which are re- 
distributed within the surrounding gas by means of a diffu- 
sive process, similar to the diffusive approximation adopted 
for the thermal conduction. In practice, the amounts of 



heavy elements (namely, Fe, Mg, O, Si, C, and global metal- 
licity Z) ejected by each SSP, i.e. star particle, is com- 
puted at each t ime-step by means of the P adova evolutionary 
tracks and the lGreggio fe Renzinil (1 19831 ) rates based on the 
number of SN events and the amount of metals released by 
the "dying" SSPs during the current time-step. The global 
chemical yield from a single SSP during each time-step is 
given by the elements released by dying (exploding) stars 
during that time-step. These are the sum of three contri- 
butions: (i) the elements locked in the dying stars at the 
moment of their birth, inherited from the proto-stellar gas; 
(ii) the heavy elements created in the core by thermonuclear 
reactions and ejected in Supernova II explosions; (iii) the el- 
ements created in the core of binary stars which end their 
life as Supernovae la. The reader should refer to iLia et all 
(2002) for a detailed description of the method. 

Numerically, the heavy elements released by a SSP are 
assigned to the gaseous mass ejected by the same SSP during 
the current time-step and subsequently spread over neigh- 
boring particles by means of a diffusion algorithm based on 
the Fick's law in spherical symmetry: 

The formalisms is the same we have adopted for ther- 
mal conduction due to diffusion of electrons. The coefficient 
Ddiff is obtained from the diffusion velocity of electrons by 
iMonaghanl (|l992r i. 



2.8 Re-ionization and the first generation of stars 

Cosmic re-ionization, as inferred from the WMAP-5 data, 
is currently believed to have taken place at z < 10. In our 
models, we include a crude treatm ent of photo-heating from 
re-ionized free electrons. Following IPawlik fe Schava (2009), 
we increase the temperature of all gas particles fueling their 
internal energy reservoir with an additional budget of 2 eV 
per proton, which is instantaneously injected at z = 9. The 
effect of this energy boost is to raise the temperature of the 
low density particles to some ~ 10 4 K. In high density re- 
gions, however, radiative cooling is so efficient to keep the 
temperature well below this value. Fig.[3]compares two snap- 
shots taken from models in which this re-ionization effect is 
(is not) included. However, it is worth noting that this way 
of proceeding lets the gas particles free to cool down soon 
after the sudden injection of energy; consequently, the cos- 
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mic web may happen to get too cold at later times (this can 
be seen in Fig.lTSIof Sect. Ell- 
in this context, we also include in a very simple fash- 
ion an early population of metal-free st ars otherwise known 
as Population HI ( Pop III) stars (e.g.. iBromm et afll2002l ; 
lYoshida et al.ll2~003l and references therein for more details). 
We assume that when a metal-free gas particle is prone to 
form stars for the first time (see Sec. I2.4|l it generates only 
Pop III objects. Furthermore, this SSP is supposed to con- 
tain only massive stars, because according to current believe 
the underlying initi al mass function is heavily skewed to- 
wards massive stars l|Bromm et al.|[2~002r ). Consequently the 
Pop III generation of stars is short lived and within a few 
million years (almost immediately on galactic time scales) 
injects lots of energy into the interstellar medium by stel- 
lar winds and supernova explosions. Therefore, thanks to 
this we can simplify the occurrence of this first generation 
as follows: gas particles undergoing star formation are not 
turned into stars, but simply raised to high temperature, 
T ~ 2 x 10 7 K, and enriched in metals, Z ~ 1CT 4 . By doing 
this we do not care to track the fate of Pop III stars as far 
as their remnants (ultra-compact objects and Black Holes) 
are concerned. Since these remnants could be the seeds to 
form to a central super-massive black hole, most likely fuel- 
ing a period of AGN activity, a complete description of the 
remnants from Pop III stars is in progress. 



2.9 General remarks 

Before going into the results of our study, we shortly discuss 
here a few possible drawbacks of our models. 

First, as already mentioned, we adopt the stellar effi- 
ciency esF = 1, in contrast to what observational data and 
simple theoretical considerations would suggest. This has 
proven to cause non-negligible effects on the dynamical evo- 
lution of the systems; they will be discussed in the Sections 
below. 

Second, no treatment of the multi-phase nature of the 
interstellar medium has been included. In contrast, differ- 
ent phases naturally develops within the gas (see Sect. 13.1)) . 
Furthermore, it is long known that numerical effects like 
overcooling may arise if the gas is modeled acc ording to the 
standard SPH method (iMerlin fc Chiosill2007ft . Given these 
considerations, a more detailed treatment of the interstellar 
medium not only is urgent but may also cure parte of the 
problems encountered with overcooling. 

Third, we have not considered the presence of nuclear 
activity, neglecting both the dynamical action of a central 
super-massive black hole, and the energy feedback from an 
AGN, both of which are currently believed to play important 
roles in the evolution of a proto-galaxy. 

Fourth, we have not included other sources of feedback 
apart from stellar - i.e., magnetic fields and their pressure 
or cosmic rays. 

Finally, our cooling functions do not take into account 
the state of the Hydrogen gas; in particular, the variation 
of the molecular fraction of Hydrogen, which strongly in- 
fluences the cooling efficiency below 10 4 K in primordial, 
non-enriched gas (see Sect. 12.3)) . 



3 MODEL RESULTS 

All the simulations have been run on the HP Blade Server 
Monster with 14 hosts and 116 CPUs at the Astronomy 
Department of Padova. Each major galaxy model made use 
of four parallel CPUs, whereas companion test models of 
small size were run on single processors. On the whole, 
the simulations took ~ 100, 000 CPU hours, requiring from 
30,000 to 80,000 time-steps. There are two groups of models: 
First group: testing the models. The first group contains 
the test computations performed to assess the model re- 
sponse to varying some important parameters. They are 
calculated with a smaller number of particles (usually from 
5000 to 10000). To save computational time, not all the mod- 
els have been calculated up to the present age. However all 
of them are carried out to the stage in which the mass as- 
sembly is completed (say at least a few Gyrs) . The following 
questions are addressed: 

(i) Star formations efficiency. These are the models al- 
ready described in sect. 12.41 They are calculated to check 
the sensitivity and overall behaviour at varying the SF di- 
mensionless efficiency csf of the star formation rate. 

(ii) Stellar feedback. A model is run with no feedback 
from stars. As expected, the star formation rate grows dra- 
matically at early times (see Fig. [4] the run had to be 
stopped very soon because the continuous collapse of cold 
gas, not contrasted by the action of heating due to feed- 
back, requires smaller and smaller time-steps and therefore 
high computational costs). Our prescription for the stellar 
feedback can therefore be considered as reasonably efficient. 

Second Group: the reference models. The second group 
contains the reference set of models, with higher resolution. 
All of them are calculated well beyond the completion of 
the mass assembly, a few even down to the present age or 
very close to it. All important data referring to the models at 
the same reference redshift z = 1 are summarized in Table [2J 
whereas in Tableware displayed similar data but for the last 
computed stage of each model (different age and redshift). 

Results for the reference models. In the following we 
shortly describe the main results for the twelve reference 
models. We start showing in Fig. [5] a snapshot of their ge- 
ometrical structure projected on the XY plane at the last 
computed. In most cases the spatial distribution of stars 
resembles that of classical ETGs; in one case a filamen- 
tary substructure, and perhaps satellites companions are 
visible. Fig. [5] displays the star formation histories of the 
twelve models along their lifetime, down to the l ast simu- 
lated epoch. As found bv lChiosi fc Carrard (|2002T ) and ex- 
pected here, the kind of star formation at work changes with 
the galaxy mass and initial density. It is monolithic-like (a 
single dominant initial episode) in high mass and/or high 
initial density proto-galaxies, it gradually turns into a broad 
rather long and mild activity at decreasing mass/and/or ini- 
tial density, and finally it gets like a series of discontinuous 
very episodes (bursting mode) for low mass and/or low ini- 
tial densities galaxies. The point will be developed in more 
detail below. 



3.1 Dynamical evolution and final morphologies 

All the models go through an initial phase of expansion, as 
they follow the Hubble flow. The central regions, in which 
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Figure 5. The spatial distribution of the stellar content of the model galaxies at the redshift of the last computed stage, projected on the 
XY plane. Left to right, top to bottom: HDHM (z = 0.22), IDHM (z = 0.77), LDHM (z = 0.49), VLDHM (z = 0.83); HDIM (z = 1.0), 
IDIM (z = 0.75), LDIM (z = 0.58), VLDIM (z = 0.15); HDLM (z = 0.36), IDLM (z = 0.22), LDLM (z = 0.05), VLDLM (z = 0.0). The 
coordinates X and Y are in proper Mpc. 



the density peak is more pronounced, soon detach from the 
outward motion and start re-collapsing, forming the core of 
the virialized structure. As time goes on, more and more ex- 
ternal regions turn around and fall onto the central virialized 
core. The process can be roughly described as a secondary 
infall, with time-scales strongly dependent on the proper- 
ties of the particular halo under consideration. More massive 
haloes globally virialize later with respect to smaller ones, 
as expected in the hierarchical cosmology. However, denser 
haloes also re-collapse sooner than the shallower ones, so the 
epoch of virialization is decided by the interplay between 
these two effects. At times depending on the depth of the 
central potential well, the central regions of the haloes be- 
gin to reach the critical density p, at which the gas becomes 
able to form stars. Due to the strongly peaked density per- 
turbations, all models soon develop a central stellar system; 
however, in the very first stages of the evolution we note the 
formation of many small clumps of stars, which soon merge 
into a single entity. In most cases, additional smaller stellar 
systems form in regions far from the center because of local 
concentrations of cold gas, eventually merging onto the ma- 
jor body of the proto-galaxy. In model galaxies of low mass, 
sparse isolated stars form within clumps of expanding gas 
(wind), forming extended haloes of extremely low density 
around the central body of the proto-galaxy. 

The whole history of the building up of the stellar com- 
ponent of the model galaxies is shown in series of panels 
in Fig. [7] Each panel show the projection on the XY plane 



(coordinates in proper Mpc) of the position of the star par- 
ticles. Each row corresponds to a total mass (high, interme- 
diate and low starting from) and each columns to an initial 
over-density: from high (left) to ver low (right). Finally each 
group contains four subpanels depending on the redshift at 
which the snaptshot are taken: from left to right and top to 
bottom the redshift is z=11.5 (7 for the very low density, 
high mass case), 5, 3, 1. 

All models finally develop a central dense stellar sys- 
tem, located at the centre of the collapsed halo of DM. Some 
amount of spurious global bulk motion of the system (due to 
numerical instabilities) can be observed in some cases (see 
Figs. [5] and [7J. At z = 1, most of the systems have almost 
completely relaxed into stable configurations, although in 
the massive ones pronounced irregular features are clearly 
visible (arms of star forming regions, diffuse haloes of stars 
around the central object). These are not formed by old pop- 
ulations infalling onto the central object, but by young SSPs 
escaping from the galaxy (this can be inferred by checking 
the birth age of the star particles and their velocities). The 
bright, spherical stellar clump visible in the bottom of the 
VLDHM model in Fig. [5] has mass M* ~ 1O 1O M , much 
larger than the masses of the low mass models. Such stellar 
clumps must therefore be considered tidal galaxies, formed 
by gas condensation. Many examples of such systems are 
known to exist (see e.g. lElmegreen et al.ll2007h . even if they 
are generally thought to be generated within interacting sys- 
tems. In our cases, local condensations of gas, compressed 
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Figure 6. Star formation histories for the twelve reference models. Left to right: high density, intermediate density, low density, very 
low density. Top row: high mass; intermediate row: intermediate mass; bottom row: low mass. 



and pushed outwards by the pressure exerted by the inner 
regions of the system, collapse to form stars while being 
ejected from the galaxy. The tidal objects are finally en- 
gulfed by the massive galaxy, but usually leave a faint ex- 
panding trail of debris. 

On the other hand, low mass systems are surrounded by 
extended haloes of sparse stellar particles, as shown in Fig. 
[5] for model HDLM at z = 1. Perhaps, these may be identi- 
fied w ith the observed intra-cluster light (ICL. lMonaco et alj 
2006), which could therefore be the result of diffuse star for- 
mation within gas escaping from low mass galaxies, rather 
than the consequence of tidal stripping and/or disruption of 
bounded systems. 

It is worth recalling here that a different choice for the 
dimensionless efficiency esF may have lead to different dy- 
namical histories, as discussed below. In particular, a lower 
efficiency could have reduced the formation of diffuse stellar 
haloes. 

All galaxies have triaxial shapes; their axis ratio at 
2 = 1 and at the final redshifts, obtained from their pro- 
jection on the XY plane, are summarized in Tables [2] and 
[21 If observed from this angle of view, all the models would 
be classified as E3 - E5 galaxies in the Hubble system. Ob- 
viously, more spherical shapes could be obtained changing 
the view angle, and almost all models could be morphologi- 
cally classified as almost perfect spherical galaxies choosing 
a suitable projection. Noticeably, the minimal initial rigid 
rotation can hardly be considered the main cause for the 
pronounced triaxiality of the models, which must be there- 
fore a consequence of their dynamical history. In particular, 



observing the first stages of the formation of these systems 
it appears that their collapsing motions took place along 
the filamentary structures of the cosmological over-densities, 
and stars tend to conserve elongated orbits. This view is 
strengthened by the fact that reducing e$F, by requiring 
a higher degree of dissipation and condensation of gaseous 
particles, yields more concentrated and spherically-shaped 
structures. It would be worth investigating how this would 
affect to present day structure these systems. Perhaps, in 
this context, a non- negligible role could also be played by the 
absence of large scale tidal motions caused by the adopted 
prescription for the evolution of the whole proto-galaxy (the 
initial sphere) in void. 

The injection of energy from stars into the ISM causes 
turbulence and large scale motions in the gas (see for in- 
stance the situation illustrated in Fig. [2]). Depending on the 
particular circumstances under which it takes place, the stel- 
lar energy feedback can lead to the ejection of large amounts 
of gaseous material from the potential well of the galaxy, or 
to the formation of condensed structures, resulting in new 
star forming sites. The minor impact of feedback in low mass 
haloes ultimately results in the large differences in their star 
formation histories, as explained below. 

Different phases naturally develop in the gas during the 
first stages of the galaxy formation as clearly shown the 
phase diagram (temperature vs. density of the gas) displayed 
in Fig. [9] for the prototype model IDIM. First, the initially 
cold, expanding gas begins to heat to a typical tempera- 
ture of T~ 10 4 K in the central regions, as soon as they 
detach from the Hubble flow and collapse. When the radia- 
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Figure 7. Evolution as a function of the redshift of the twelve models with different total mass and initial densities. The positions of 
the star particles are projected on the XY plane (proper Mpc). The upper row corresponds to the high mass case, the central row to 
the intermediate mass case, and the bottom row to the low mass case. The initial density decreases from left to right. So the models 
displayed in the three groups of small panels in the left column are HDHM, HDIM, HDLM. The same in the other columns but for 
the ID-, LD- and VLD- cases, respectively. In each sub-panel the redshift changes from left to right and top to bottom: z ~ 11.5 (7 for 
VLDHM), 5, 3, 1. 



tive cooling becomes effective, this gas begins to cool and 
to form stars. The ejecta from the SSP subsequently heat 
the surrounding^. Now, the infalling gas is forced to inter- 
act with these heated regions. If the hot gas is enough, this 
may lead to expanding shock fronts (see below). Thus, after 
several Myr of activity, the gas can be roughly split as fol- 
lows. A cold and dense phase is formed by gas in the central 
regions of the galaxies, where stars are forming. On the op- 
posite extreme, hot and thin gas either is leaving the halo 
(galactic wind) after being heated by stellar feedback, or it 
is infalling onto the central object, after being shock-heated 
because of the interaction with the central, hot regions. Two 
more phases are clearly detectable: the hot and dense gas, 
where the stellar feedback is heating the gas at very high 
temperatures before these regions expand under the action 
of the local pressure, and the cold but thin gas belonging 
to the cosmological cosmic web (before re-ionization takes 
place, or cooling after re-ionization; see Sect. 
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Figure 8. Low mass models develop extended haloes of stellar 
particles around the central object (in this picture, model HDLM 
is shown at z = 1; distances are in proper Mpc). 



3.2 Star formation histories 

Looking at Fig. [6] it is clear that the initial conditions of 
each halo deeply influence the star formation history (SFH) 

3 Once more, it is worth pointing out that heating of the sur- 
rounding gas is due to the mechanical pressure exerted by the 
stellar particles, and not to a redistribution of their thermal bud- 
get; see Sect. 12.51 



of the parent proto-galaxy. Even at a first glance, two clear 
trends can be easily identified. 

First: the more massive the halo, the more peaked is the 
SFH. All massive haloes M to t — 10 13 M Q ) show a single, 
strong burst at very early times, followed by a slow decrease 
down to quiescence. A weak tail of activity persists down to 
later times in the less dense systems, thus accounting for a 
weak dependence on its density (or environment). Smaller 
haloes, in turn, display a more fragmented and prolonged 
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Table 2. Properties of the twelve major models at z = 1 (i.e. after ~ 6 Gyr of evolution). Left to right: total stellar mass, total mass 
of the virialized halo, fraction of the total initial mass composed by star, fraction of the initial gas converted into stars, virial radius 
of the whole system (i.e. stars, Dark Matter and gravitationally bound gas; it is the radius at which the mass density is 200 times the 
background density at the epoch of virialization, which we assume to have happened at t ~ 2 Gyr looking at the energy trends of the 
models: see Fig. 1221 1. half mass radius of the stellar system projected on the XY plane, axis ratios of the stellar system projected on the 
XY plane. 



Model 


Ah [M ] 


M vir [Mq] 


M,/M vir 


Af* / M g as,ini 


r v ir,tot [kpc] 


n/2,* [kpc] 


b/axY 


HDHM 


7.8 x 10 11 


1.1 x 10" 


0.071 


0.27 


153.0 


15.7 


0.53 


IDHM 


7.5 x 10 11 


1.0 x 10" 


0.075 


0.26 


141.8 


16.8 


0.54 


LDHM 


7.4 x 10 11 


9.4 x 10 12 


0.080 


0.26 


133.8 


14.5 


0.57 


VLDHM 


6.3 x 10 11 


7.5 x 10 12 


0.080 


0.22 


112.5 


11.2 


0.52 


HDIM 


2.0 x 10 10 


1.6 x 10 11 


0.13 


0.45 


37.6 


5.7 


0.62 


IDIM 


1.9 x 10 1U 


1.5 x 10 11 


0.13 


0.43 


35.7 


5.8 


0.63 


LDIM 


1.9 x 10 1U 


1.4 x 10 11 


0.14 


0.42 


33.3 


5.3 


0.58 


VLDIM 


1.3 x 10 1U 


1.1 x 10 11 


0.12 


0.29 


28.3 


5.8 


0.57 


HDLM 


1.2 x 10 s 


2.6 x 10 y 


0.028 


0.17 


9.2 


2.2 


0.69 


IDLM 


1.0 x 10 s 


2.5 x 10 a 


0.040 


0.14 


10.0 


2.5 


0.64 


LDLM 


8.9 x 10 7 


2.3 x 10 a 


0.021 


0.13 


11.8 


2.2 


0.56 


VLDLM 


5.0 x 10 7 


1.7 x 10 y 


0.011 


0.07 


10.5 


2.7 


0.58 



Table 3. Properties of the twelve major models at the final stage of the simulated evolution. Left to right: final redshift, final age, total 
stellar mass, total mass of the virialized halo, fraction of the total mass composed by star, fraction of the initial gas converted into stars, 
virial radius of the whole system (i.e. stars, Dark Matter and gravitationally bound gas; it is the radius at which the mass density is 200 
times the background density at the epoch of virialization, which we assume to have happened at t ~ 2 Gyr looking at the energy trends 
of the models: see Fig. |22> . half mass radius of the stellar system projected on the XY plane, axis ratios of the stellar system projected 
on the XY plane. 



Model 




tiast [Gyr] 


M, [M Q ] 


M vir [Mq] 


M,/M vir 


JVf* / Mgas,ini 


r v ir,tot [kpc] 


r i/2,* [kpc] 


b/axY 


HDHM 


0.22 


11.0 


7.5 x 10 11 


1.5 x 10 13 


0.050 


0.26 


153.0 


15.6 


0.56 


IDHM 


0.77 


8.0 


7.4 x 10 11 


1.5 x 10 13 


0.050 


0.26 


141.8 


16.5 


0.48 


LDHM 


0.50 


8.7 


7.3 x 10 11 


1.5 x 10 13 


0.049 


0.25 


133.8 


15.8 


0.57 


VLDHM 


0.83 


6.6 


6.3 x 10 11 


1.3 x 10 13 


0.048 


0.22 


112.5 


11.2 


0.52 


HDIM 


1.0 


5.8 


2.0 x 10 iu 


2.1 x 10 11 


0.10 


0.45 


37.6 


5.7 


0.62 


IDIM 


0.75 


7.0 


1.9 x 10 1U 


2.1 x 10 11 


0.08 


0.43 


35.7 


5.8 


0.63 


LDIM 


0.58 


8.1 


1.9 x 10 1U 


2.0 x 10 11 


0.10 


0.42 


33.3 


5.2 


0.75 


VLDIM 


0.15 


11.8 


1.7 x 10 10 


1.4 x 10 11 


0.12 


0.38 


28.3 


4.9 


0.83 


HDLM 


0.36 


9.7 


1.5 x 10 s 


3.3 x 10 a 


0.045 


0.19 


9.2 


2.3 


0.74 


IDLM 


0.22 


11.0 


1.4 x 10 s 


3.3 x 10 y 


0.04 


0.16 


10.0 


2.4 


0.67 


LDLM 


0.05 


13.0 


1.4 x 10 s 


3.2 x 10 y 


0.04 


0.19 


11.8 


2.1 


0.79 


VLDLM 


0.0 


13.7 


1.0 x 10 s 


3.0 x 10 y 


0.03 


0.10 


10.5 


2.7 


0.65 



Table 4. Properties of the star formation history for the twelve models. Left to right: maximum peak of star formation rate [M0/yr], 
epoch of maximum activity [Gyr], percentage p* of the final stellar mass assembled at four different redshifts, and epoch at which different 
fractions / of the final stellar mass (i.e. at z = 1) has formed [Gyr] (with corresponding redshift). 



Model 






P*,z=10 


P*,z=5 


P*,z=2 


P*,z=1.5 


%=l/2 


(t, z)/=99/100 


HDHM 


1380 


0.5 


1 4.9 


67.3 


99.5 


100.0 


0.863 6.5 


3.029 2.2 


IDHM 


1020 


0.6 


2.9 


55.2 


98.0 


100.0 


1.097 5.4 


3.778 1.8 


LDHM 


800 


0.9 


<0.1 


39.7 


93.2 


98.5 


1.417 4.4 


4.398 1.5 


VLDHM 


350 


1.5 


0.0 


5.2 


67.2 


80.1 


2.423 2.7 


5.083 1.2 


HDIM 


19.0 


1.2 


2.8 


35.5 


92.2 


98.9 


1.589 4.0 


4.392 1.5 


IDIM 


13.5 


1.5 


0.1 


19.9 


86.0 


98.3 


1.976 3.3 


4.530 1.4 


LDIM 


12.0 


2.5 


<0.1 


9.5 


69.4 


85.9 


2.568 2.6 


5.806 1.0 


VLDIM 


8.0 


2.5 


0.0 


1.1 


42.8 


66.3 


3.537 1.9 


5.837 1.0 


HDLM 


0.135 


1.5 


2.3 


17.9 


61.5 


76.1 


2.780 2.4 


5.940 1.0 


IDLM 


0.075 


3.0 


0.6 


10.2 


49.1 


69.6 


3.480 1.9 


5.850 1.0 


LDLM 


0.06 


2.0 


0.2 


6.2 


54.9 


79.2 


3.000 2.2 


5.800 1.0 


VLDLM 


0.55 


3.8 


0.0 


4.8 


25.2 


46.0 


4.400 1.5 


5.800 1.0 



star formation activity, often alternated with periods of qui- 
escence. 

Second: lower densities correspond to delayed activity 
(the first burst starts later) and minor overall efficiency (the 
peaks of activity have lower magnitude). A delay up to ~ 1 



Gyr in the beginning of the SF activity can be seen in the low 
density models. Moreover, the maximum rate is generally 
halved, and the activity is prolonged (this second effect is 
of crucial importance in intermediate mass haloes, M to t — 
1O 11 M0). Apparently, the mass of the halo is the dominant 
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Figure 9. Phase diagram density p (g/cm 3 ) vs. temperature T 
(K) of gas particles for the IDIM model, at four stages of evolution 
(z ~10.0, 7.0, 5.0, 3.5). The hot gas decreasing its density is 
leaving the galaxy, due to stellar feedback. 
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Figure 10. Star formation histories of the whole galaxies (solid 
black line, large circles) compared with the corresponding ones 
within the central region (ri nner ~ 20 kpc; red dashed lines, large 
points) for two large mass systems. Left panel: model HDHM. 
Right panel: model VLDHM. 



characteristic to determine its evolution, but intermediate 
mass haloes can have very different histories depending on 
their global over-density. 

Table [4] summarizes the properties of the SFHs. In- 
teresting quantities are the intensity of star formation (in 
Mq/yt) at its peak value and the time (in Gyr). Further- 
more, the gradual building up of the stellar mass in a galaxy 
is measured by the percentage p* of the total stellar mass 
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Figure 11. Same as in Fig. 1101 but for two small systems (here 
dinner — 2 kpc). Top panel: model HDLM, bottom panel: model 
VLDLM. 



at redshift z — 1 assembled at selected values of the redshift 
(z=10, 5, 2 and 1.5), and the time (redshift) at which 50% 
and 99% of the stellar mass at z=l is assembled. This is 
indicated by the pairs (t, z) /=i/2 and (t, z) / = gg/ioo in Table 

a 

All the models begin forming Population II stars be- 
tween z ~ 20 and z ~ 7. The HDHM halo experiences 
the strongest and earliest peak of activity, reaching some 
1.5 x 10 3 M©/yr at z ~ 10. More than 50% of its total stel- 
lar mass is assembled before z ~ 6. If observed now from the 
Earth, this galaxy would appear massive (M* — 3x 10 M©) 
and essentially red at z ~ 4. Its tail of SF activity contin- 
ues down to 2 ~ 1.6, contributing 19% of the total stel- 
lar mass. Looking at less dense haloes of the same mass, it 
can be noticed how the maximum of activity decreases in 
magnitude and is delayed to later times, down to a peak of 
7.5 x 1O 2 M /yr at z ~ 5 in the VLDHM model, where more- 
over an intense activity is prolonged to much later epochs. 
Despite the similarity between the final objects formed from 
massive haloes, a clear trend can be extrapolated from the 
data given in Table [4] the lesser dense the halo, the later 
the mass assembly takes place (with the VLDHM halo as- 
sembles half of its total stellar mass at z ~ 2.7). The overall 
efficiency of the star formation process in these models is 
20%-30%. 

Intermediate mass systems apparently have more com- 
plex histories. Moving from denser to less dense systems, 
the shape of the star formation curve changes dramati- 
cally. Strongly over-dense haloes undergo an intense burst 
at early times, much as the most massive systems. However, 
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Figure 12. Date of birth of the stars (in Gyr) as a function of 
their radial positions (distance from the barycenter in kpc) for 
the high and low mass models models (left to right: decreasing 
density; top to bottom: decreasing mass), at the final time-step 
of their evolution (the same age as in Fig. [5J . The radial distance 
is calculated from the barycenter of the stellar system (kpc). The 
virial radii of the stellar systems are similar to the extension along 
the X axis. 



less dense haloes gradually display more prolonged, and un- 
even histories, with many different peaks. Apparently, in this 
range of masses the initial density plays a dominant role in 
deciding the fate of the mass assembly of the galaxy. The 
efficiency of the star formation process is very high in these 
cases: 30% -50% of the initial gas is converted into stars. 

All low mass systems show a continue, bursty SFH, 
sometimes with long periods of quiescence. In these systems, 
the feedback from stars is sufficient to halt the star forma- 
tion process soon after its beginning. However, this means 
that as soon as the first stars cease their strong feedback 
activity, the gas is able to re-collapse and form new genera- 
tions of stars. This galactic "breathing" leads to systems in 
which many different generations of stars are present. The 
overall efficiency of the star formation mechanism is lower 
in this cases, with 10-20% of the initial gaseous mass turned 
into stars (giving high M / L ratios, comparable to those of 
many observed dwarf galaxies). 

Figs. [10] and [TT] compare the SFH of four models 
(HDHM, VLDHM, HDLM, VLDLM) across the whole sys- 
tems with the SFH within their central regions (that is, 
within the central 20 kpc in massive systems and within 
the central 2 kpc in the small ones). It can be noticed that 
massive systems build their central regions before their out- 
skirts, while in less massive galaxies the SF process is active 
both in the inner and in the outer regions during the whole 
lifetime of the galaxies . This is consistent with some recent 
observational findings dvan Dokkum etaitold) on the way 
galaxies are built-up. 

3.3 Stellar populations: ages, mean metallicities, 
and gradients in metal content 

Stellar ages. The distribution of the galactic ages (Tg,») at 
which the star particles are created as a function of the ra- 
dial distance is shown in Fig. 1 121 limited to the high and low 
mass models (the behaviour of the intermediate mass ones 



0.2 


-0.2 

a" -o.4 

o 

-0.6 
-0.8 




8 9 10 11 12 

log[M*] 

Figure 13. Averaged metallicity log[Z/Zo] versus stellar mass 
logAf* [Mq] for the twelve reference models, compared to 
the observed mass- metallicity relation. Black lines: data from 
iGallazzi et al. (2005) (solid line: median distribution of their sam- 
ple; dashed lines: 16th and 84th percentiles). Red open circles: 
averaged metallicity of all stellar particles in each model. Blue 
open circles: averaged metallicity within the inner 5 kpc in each 
model. 
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Figure 14. Gradients in metallicity. Left to right, top to bot- 
tom: models HDHM, LDHM, HDLM, LDLM at their last com- 
puted snapshot. Plotted is the metallicity log[Z/ZQ] versus the 
fractional radial distance r/r e of all particles. Also shown are the 
median-binned metallicity profile, the number of particles being 
the same in each bin (blue and red solid lines, referring to the in- 
ner and outer region, respectively), and the linear best fit (black 
solid line; the slope is written in the left bottom corner of each 
panel). See text for details. 



falls in between these two). The real star particles ages are 
T» = Tg — Tc,i- In early epochs, the stars are preferentially 
created in the central regions, then the star forming activity 
expands to larger radii (inside-out mechanism), and mov- 
ing towards the present time, the stellar activity tends to 
shrinks again towards the centre. This simply mirrors the 
star formation history and the mechanism of mass assembly 
presented above. 

Mean metallicity. Turning now to the metal content 
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reached by the model galaxies, we look at the average stellar 
metallicity as a function of the total stellar mass. The re- 
sults are shown in Fig.[T3]where the red circles refer to all the 
star particles, whereas the blue ones only to those within the 
central 5 kpc of each model. The black lines are the median 
(solid) and the 16 th and 84th percent iles (dashed) of the 
data presented in (jGallazzi et al.ll2005l ). There is marginal 
agreement between data and theoretical results. The abso- 
lute values of the theoretical metallicities are too low com- 
pared with the observational estimates, in particular for the 
most massive models which are almost an order of magni- 
tude too metal-poor. This may be due to several factors. 
First, the theoretical values are straightforwardly computed 
averaging the metallicity of all particles within the chosen 
radius without any further refinement; of course, observed 
metallicity are luminosity weighted, while in our models the 
contributions from very metal-poor, primordial (old) and 
faint stars the same weight of more recent stellar popula- 
tions. Second, part of the discrepancy may arise because of 
the different methods used to obtain observational and the- 
oretical metallicities. For example, the use of the Kroupa 
IMF, which is tailored to fit the solar vicinity, may lead to 
underestimate the enrichment by massive stars in other en- 
vironments. Recent studies seem to indicate that the IMF 
in massive ETGs is more skwed toward m assive stars than 
commonly assumed (|Cappellari et al1l2012l ) thus implying a 
net increase of metal production by massive stars in these 
systems. Third, the high star formation efficiency can again 
play a role, favouring the formation of low metallicity stars 
at early times and reducin g the average me tallicity of the 
m odels. Finally, as shown bv lScodeggiol l|200ll ) and discussed 
bv lTremonti et al. (2004), at least some part of the observed 
color- (and hence metallicity-) mass relation might just be a 
spurious consequence of an aperture effect, as observed col- 
ors (and metallicities) are usually measured within a fixed 
aperture, rather than a given fraction of the galaxy effective 
radius, so that the absolute values are not well determined. 
In conclusion, considering that chemical enrichment in NB- 
TSPH simulation is still far from being fully satisfactory, we 
are inclined to say that theory and data marginally agree 
each other. 

In passing, we point out that the models naturally 
reproduce the observed trend of the a-enhancement phe- 
nomenon in massive galaxies. The early peak of star forma- 
tion in high mass systems, followed by quiescence, ensure 
the enhancement in light elements, whereas low mass sys- 
tems with prolonged SFH are less enhanced, as it should be 
(see the discussion in iTantalo fc Chiosil200a) . 

Metallicity gradients. Radial gradients in spectro- 
pho tometric properties of E TGs are known to exist (see 
e.g. lLa Barbera et all |2010| . and references therein). The 
majority of ETGs in the local U niverse features red cores 
l|La Barbera fc de Carvalhol 120091). and less freq uently blue 
cores, mostly found at low mass liuh et al.l201Ch . The radial 
gradients are interpreted as age and me tallicity gradien ts of 
the constituent stellar populations (e.g. Gonzalez 1993). 

In Fig. [TJ we plot the metallicity of each particle as a 
function of projected distance to the galaxy center (from left 
to right, and top to bottom, models LDHM, HDHM, LDLM 
and HDLM are displayed). For each model, the median- 
binned metallicity profile is shown, the number of particles 
being the same in each bin. In order to estimate the metallic- 



ity gradient, we performed an orthogonal least-square fit of 
the metallicity profile in the range 0.1 r e to 1 r e , i.e. the sam e 
range used in observational studies (|Peletier et alj Il990h . 
The fitting was repeated for 100 different two-dimensional 
projections, the average slope of the best-fitting lines giving 
the metallicity gradient, Vz = dlog Z/dlog R, of a given 
model. We repeated the computation of Vz by applying dif- 
ferent weighting schemes, where each particle was assigned 
no- weight, mass- weighted, and luminosity- weighted, finding 
no significant differences within the errors, especially at high 
mass (see Tab. [5}. 

Looking at the results in some detail, high-(low-) mass 
models exhibit metallicity gradients spanning the range from 
-0.4 (-0.02) to -0.27 (0.16), depending on the weight- 
ing scheme and environment. In general, low-density models 
tend to have more negative metallicity gradients. Moreover, 
at given density, high- (relative to low-) mass systems have 
more negative gradients. 

In general, the majority of the stars in the central re- 
gions of the massive systems have higher metallicity than 
those in the outskirts; on the contrary, in less massive sys- 
tems the metallicity of the central stars is generally equal to 
that of the peripheral ones. Even if they are not straightfor- 
wardly visible in Fig. 1141 interesting features to recall are: 
(i) on the average, the metallicity distribution of the star 
particles gets narrower at increasing radial distance; (ii) the 
maximum values are always reached in the centre, (iii) great 
deals of stars are made of recycled gas; (iv) finally the sim- 
ulations often show the presence tidal satellites which have 
their own chemical history. 

Although the metallicities of our models do not match 
observations in a absolute sense, the relative variations of 
metallicity (between, e.g., different masses, and/or different 
galacto-centric distances) are likely more robust, as they 
are less dependent on the absolute calibration of the in- 
gredients the models rely on. This motivates for a com- 
parison of the predicted Vz's with observational results. 
A steepening of the metallicity gradient with galaxy mass, 
qualitatively consi stent with that seen in our models, has 
been reported by lLa Barbera et al] (l2010h. and, i n depen - 
dent of galaxy environment, by lLa Barbera et al.l <|201ll ). 
On the other hand, other observational studies, based on 
smalle r samples, do not fi nd a clear correlation of Vz with 
mass aor et aill2009l ). or even a double- value behavior 
l|Tortora et alj201Ch . For a galaxy mass of 3 x 10 11 M© (sim- 
ilar to that of ~ 7 x 10 11 M Q of models HDHM and LDHM), 
La Barbera and collaborators found metallicity gradients of 
-0.37 ± 0.02 and -0.41 ± 0.02 for low- and high-density 
ETGs, respectively. At low density, the estimate of is fully 
consistent with our model predictions (for no-weight and 
mass- weighted Vz's). At high-(relative to low-) density, the 
models exhibit shallower gradients, by ~ 0.07, in contrast 
to the data, where a (marginal) difference of —0.04 ± 0.03 
is found. For what concerns low-mass models, they exhibit 
almos t null gradient s , cons istent with the observational find- 
ing of lSpolaor et~ail (|2009l ) that the Vz vanishes at a galaxy 
mass of ~ 10 9 M© . Whilst these results point to a fair agree- 
ment between our models and (some) available data-set, we 
postpone a more thorough investigation of both metallicity, 
age, and color gradients, as well as their correlation with 
other model properties, to a forthcoming contribution. 
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Table 5. Metallicity gradients of four ETG models in the radial range of 0.1 r e to lr e . Column 1 is the model ID label. Columns 2, 
3, and 4, correspond to cases where the Vz is computed with no-weight, mass-weight, and luminosity-weight assigned to each particle. 
Errors are the rms of Vz estimates among 100 projections of each model. 



ID 


Vz (no-weight) 


Vz (mass-weighted) 


Vz (lum-weighted) 


(1) 


(2) 


(3) 


(4) 


LDHM 


-0.39 ± 0.04 


-0.40 ±0.04 


-0.33 ± 0.03 


HDHM 


-0.32 ±0.04 


-0.31 ± 0.04 


-0.27 ±0.03 


LDLM 


0.00 ± 0.04 


+0.01 ± 0.02 


0.00 ±0.06 


HDLM 


+0.09 ±0.02 


-0.02 ± 0.01 


0.16 ±0.03 



3.4 Cold and hot gas accretion and ejection 
(galactic winds) 

Fig. 1151 shows the gas temperature vs. density diagram for 
models of our set at z = 2. As expe cted from theoretical 
considerations (jRees fc O strikcr 1973), a bimodal gas accre- 
tion can be observed in our models, as a function of the total 
mass of the proto-galactic halo. In low mass haloes, the gas 
accretes from the cosmic web without being shock heated to 
the virial temperature, and proceeds to flow along filaments 
towards the center of the halo, where it will eventually shock. 
In contrast, massive haloes soon shock the infalling gas to 
very high temperatures, suppressing the formation of local 
substructures and forming a central distribution of hot gas. 
haloes which do support shocks close to the virial radius are 
expe cted to contain a quasi-h ydrostatic atmosphere of hot 
gas. iBirnboim fc Dekell (|2003h claim that shocks can only 
form close to the virial radius in haloes with mass greater 
than 1O 11 M0 for primordial gas (or around 10 12 Mq for gas 
of Solar metallicity), in good agreement with our models. 

Looking at Fig. 1151 it can be noticed that the situation 
shows a strong dependence on the mass of the halo. In mas- 
sive and dense systems, almost no cold gas (T< 10 4 K) is 
present, so the star formation cannot proceed further. More- 
over, an expanding shock front is evident (most of the gas 
is gathered in a narrow density strip moving outwards at 
high temperatures); this will result in a hot galactic wind. 
The hot but thin gas is being shocked during its infall, in 
a hot accretion mode. The low density, cold tail is formed 
by gas expanding and never collapsing onto the galaxy, cool- 
ing down after being photo-heated by reionizing radiation at 
z = 9 (see Sect. I2.8|) . More cold gas is present in less dense 
systems, both in high density regions (potentially forming 
new stars) and in the outskirts, the shock front being less 
pronounced. A large amount of dense (p ~ 10~ 25 g/cm 3 ) 
and hot (T> 10 6 K) gas is also present in massive systems 
which are still forming stars, due to the action of feedback. 

In the intermediate mass systems, a large amount of 
cold and dense gas is still present, ensuring a continuing star 
formation. In low mass systems, almost no gas is heated to 
T> 10 5 . Thus, the gas is ejected from the galaxy at low 
temperatures, and no shock front develops. 

Figs. [16] and [17] plot the state of the gas in two ref- 
erence models, namely HDHM and LDLM; they show the 
radial velocity as a function of the radial distance at three 
different epochs and compared to the escape velocity (the 
blue solid lines; see the caption for more details). The mass 
of the ejected gas increases with time in all the models. The 
average velocity is proportional to the mass of the system, 
with massive galaxies generating winds at ~ 2000 km/s, 
and low mass galaxies at about one tenth of the speed. In- 



terestingly, the average radial velocity is proportional to the 
distance, like in a Hubble-like flow; the ejection of the gas 
is therefore an explosive phenomenon, rather than a slow 
outflow. However, looking more carefully it can be noticed 
that there are different linear velocities at the same radial 
distances, meaning that many explosive events have taken 
place. There is also a huge mass of gas moving outwards 
at slower speed, at late times (particularly in the HDHM 
model). 



In Fig. [I8]we plot fractional amounts of escaping gas as 
a function of time, for four of the models (namely, HDHM, 
LDHM, HDLM, LDLM; the other models show similar be- 
haviours). We plot the trends only after the systems have 
virialized; in other words, the total gravitational energy 
must have become constant to avoid the initial phases of 
galaxy assembly during which the gas accumulates toward 
the center and the galaxy is subject to episodes of quick ex- 
pansions and/or contractions. The red dotted lines display 
the fractional mass of gas which is found outside a given ra- 
dius r e xt from the centers of the systems at the given time. 
For massive systems we adopt r ex t = 200 kpc, while for small 
systems r ex t = 10 kpc (these values are similar to the virial 
radii of the models, see Tab. [3]). The blue dashed lines are 
relative the fraction of gas which has an outward-directed 
radial velocity, and is moving faster than the Hubble flow at 
the correspondent redshift. Finally, the black solid lines are 
the fraction of gas which is also above the escape velocity of 
the system. Clearly, two different regimes are at work. In the 
high mass realm, the fraction of gas which is really leaving 
the galaxy is only a few percent. On the contrary, in low mass 
systems more than half of the total gas has sufficient velocity 
to escape. This is easily understandable. Assuming constant 
density in models of different total masses, the gravitational 
potential energy of a system is proportional to M 5 ^ 3 , while 
the energy injection from SN explosions is proportional to 
the number of stars in the galaxy, i.e. roughly oc M. Thus, 
the binding energy of a system grows more rapidly than the 
kinetic energy of its components, making it more difficult to 
escape from the system. However, this distinction is theo- 
retical, and what we observe in the real Universe is simply 
the mass of gas that can be found outside a galaxy, moving 
outwards (i.e., the blue dashed line in the diagram). In this 
case, the situation in the two cases is much similar, however 
still with a higher efficiency in the low mass regime. The 
escaping gas is chemically enriched, as illustrated in Fig. [19] 
for the HDHM model. 
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Figure 15. Phase diagram density p (g/cm 3 ) vs. temperature T (K) of gas particles for the twelve reference models, at redshift 2 = 2. 
In each panel: left to right: high density, intermediate density, low density, very low density. Top row: high mass; intermediate row: 
intermediate mass; bottom row: low mass. 



3.5 Critical discussion of the results 

Is the strong dependence of the SFH on the initial galaxy 
mass and density a physical result, or is it just a spurious 
numerical effect? The large difference in particle masses be- 
tween the various simulations, a necessary choice to model 
galaxies spanning a wide range of sizes, is the major rea- 
son for this serious uncertainty. In high mass models, the 
gas particles have a mass vhum — 5 x 1O 7 M0, whereas in 
the low mass models the mass of a gas particle is wilm — 
1.2 x 10 4 M Q , i.e. a factor of ~ 4 x 10 3 times smaller. Does 
this difference in mass resolution affect the global results of 
our study? 

To check this issue, we run three low resolution models. 
The initial haloes for these test models were crudely ob- 
tained from the three high density models HDHM, HDIM, 
and HDLM, by simply considering one particle out four 
when reading the input files, and simultaneously multiply- 
ing the particle masses by a factor of four. Although this 
procedure is not rigorous, still it can be safely used to get a 
quick insight of the problem. The new models shortly indi- 
cated by lrHDHM, lrHDIM, lrHDLM have the total number 
of particles reduced by a factor of 4 and can be calculated 
on a single CPU. Their star formation is compared to that 
of their "parent" models and the results are plotted in Fig. 
1201 Even if some minor differences can be noticed, it is evi- 
dent that the global trends are quite similar at the different 
resolutions. This is reassuring, because it implies that we 
can trust the results of our reference models, at least as far 
as their different star formation histories is concerned. This 
also shows that meaningful results can be obtained even us- 
ing a small number of particles with much less numerical 
effort. 

Finally, we consider the effects of imposing constant 



minimum threshold values to the softening and smoothing 
lengths of particles. To do so, we re-run the low mass mod- 
els imposing minimum values e m i n = h m i n = 100 pc. We 
plot the SFHs of the resulting models in Fig. 1211 and com- 
pare them with the SFHs of the same models belong to the 
reference set in which no limitations have been imposed to 
the softening and smoothing lengths. Clearly, there are non- 
negligible differences. The global efficiency in the star forma- 
tion process is lower in the models including the threshold. 
Moreover and perhaps even more important, the dynamical 
evolution is substantially different, because in these models 
star formation is essentially halted after a strong peak of ac- 
tivity, whereas in the reference models it proceeds at nearly 
constant efficiency. To cast light on this issue we look at 
the temporal variation of the kinetic, gravitational potential, 
thermal, and total energies of the models under considera- 
tion. They are plotted in Fig. [22] It is clear that imposing the 
limit values to the softening and smoothing lengths causes 
an uncorrect behaviour at the epoch of the strong ejection 
of hot gas following the first intense burst of star formation, 
when the dynamical stability of the system is put in trou- 
bles. The model response is not physically adequate to drive 
the correct dynamical behaviour of the whole system. This 
means that the wrong choice of these parameters can deeply 
affect the overall evolution of the model. Variable softening 
and smoothing lengths are therefore highly recommended. 
Nevertheless, a great deal of the general behavior of the ref- 
erence models (in particular the gross features of the star 
formation history) still survives also in those with the con- 
strained softening and smoothing lengths, thus securing that 
even with wrong choices for these two parameters we catch 
the essence of the problem. 
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Figure 16. Galactic winds in the HDHM model. Red points: 
radial velocities of gaseous particles (km/s) against the distance 
from the barycenter of the systems (Mpc). Black circles: mean 
radial velocity in the spherical shell at radius r. Solid blue line: 
escape velocity as a function of the radius. Top to bottom: z = 4.4, 
1, 0.2. 
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Figure 18. Fractions of gas leaving the systems as a function of 
time. Top panels, left to right: HDHM, LDHM models. Bottom 
panels, left to right: HDLM, LDLM models. Red dotted lines: 
fractional mass of gas found outside r ex t from the centers of the 
systems at the given time. Blue dashed lines: fraction of gas which 
has an outward-directed radial velocity, and is moving faster than 
the Hubble flow. Black solid lines: fraction of gas which also has 
radial velocity above the escape velocity of the system. See text 
for more details. 
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Figure 17. Galactic winds in the LDLM model. Symbols as in 
Fig. 19. Top to bottom: z = 2.2, 1, 0.05. 



Figure 19. Metallicity of the ejected gas (wind) in the HDHM 
model at Zf in = 0.22. Distances are in proper Mpc. Chemically 
enriched gas (red particles), with Z ~ Zq ~ 0.019, is leaving 
the galaxy (which is located near the origin). Note the different 
metallicity, one or two orders of magnitudes lower, of the gas 
particles in the outskirts (green and blue): these do not belong 
to the "galactic wind" , since they haven't been ejected from the 
galaxy; they are part of the "cosmic web" . 



4 MASS DENSITY PROFILES 

The geometrical structure of the model galaxies is best illus- 
trated by the surface and volume density profiles. To com- 
pare these we choose a certain evolutionary stage, namely 
the model galaxies at a certain redshift, i.e. z = 1. Fig. [23] 
displays the surface density radial profiles projected on the 
XY p lane at such redshift. Over-plotted is the best ISersicl 
(1968) profile, 



w here r e is the effe ctive radius of the galaxy, as defined e.g. 
in Irlernquistl fl990), oo is the surface density at r e , and m 
is the Sersic index (if m = 4 the de Vaucouleurs profile 
is obtained). All profiles are computed starting at 0.2% of 
the virial radius of the galaxies to avoid the very central 
region of a galaxy where softening may introduce spurious 
numerical effects. The best-fitting Sersic index is m ~ 4, 
m ~ 1.5, and m ~ 2.5, for high-, intermediate-, and low- 
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Figure 20. Star formation histories for the three low resolution 
models, lrHDHM, lrHDIM, lrHDLM (in the bottom panels the 
corresponding SFH for the high resolution models are plotted for 
comparison). The rates have been normalized to the mass of the 
particles. 
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Figure 21. Star formation histories for the low mass models, 
without (top panels) and with (bottom panels) imposing a mini- 
mum softening and smoothing length. 



mass models, respectively. In other words, high- (relative to 
low- and intermediate-) mass models tend to have higher m, 
in qualitative agreement with t he existence of a luminosity- 
Sersic index relation of ETGs (|Caon et al.lll993l ). However, 
one should notice that the most mass ive ellipticals in the lo - 
cal Universe tend to have m ~ 8 (e.g. lFerrarese et al.llioOril) . 
while we find m ~ 4. Moreover, our intermediate-mass mod- 
els have somewhat lower m than the low-mass ones, which is 
only marginally consistent with luminosity-m relation, con- 
sidering its large scatter. 

The overall agreement between the models and the Ser- 
sic curves is good in the external regions. However, a clear 
departure from the expected fits is evident in all models 
at small radii (some fraction of the effective radius, indi- 
cated with the black diamond shown in each panel). The 
central regions of the model galaxies tend to flatten out with 
a plateau in the mass density profile. Given the adaptiveness 
of the force softening, this feature can hardly be ascribed to 
numerical artifacts. The close scrutiny of the problem, sug- 
gests us ascribe it to the high value of the dimensionless 
efficiency, esF ~ 1, that we have adopted for the star for- 
mation rate. 

This conclusion is strengthened by the three- 
dimensional mass density profiles shown in Fig. 1241 These 
profiles are computed starting at 0.05% of the virial radii, 
and they are compared with the best fits of the Hernquist 
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Figure 22. Energy trends in the low mass models. In each 
panel the global values of the energies (in code units) are plot- 
ted versus time, superposing the results of different models. Left 
panel: models without minimum softening/smoothing lengths; 
right panel: models with minimum softetning/smoothing lengths 
(f m ,„ = h m i n = 100 pc). Red dashed lines: kinetic energy; 
blue dashed lines: gravitational potential energy; green solid lines: 
thermal energy; black lines: total energy. 



(|Hernauistlll99b1 ) and (|Navarro et alJll996h profiles (NFW), 
for stars and DM, respectively. While DM haloes smoothly 
follow the expected NFW profiles, the density of the stel- 
lar component never exceeds that of the DM, and the cen- 
tral regions of the model galaxies are dominated by the 
DM, contrary to the expe ctations (e.g. iTortora et al. 1 120091 : 
iPadmanabhan et alj|2004 but see also Grillo, 2010). 

The relations between the efficiency of star formation 
and the gas density at which star form is straightforward. 
If star formation proceeds slowly (low values of esF, the 
gas can get high densities before it gets consumed by star 
formation itself. The opposite if star formation proceeds at 
high efficiency, high values of esF- As a consequence there 
is an immediate correlation between star formation and the 
dynamical behavior of a galaxy. Allowing gas particles to 
turn into stars as soon as they fulfill the required numerical 
criteria (that is, with unit efficiency) will generally result 
in quite diffuse stellar systems, with shallow potential wells 
and large effective radii. On the contrary, a low star forma- 
tion efficiency leads to dense clumps of cold gas which give 
birth to more compact stellar clouds (and to a delay in the 
beginning of the activity). 

Interestingly, this puts a constraint on the modality in 
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Figure 23. Stellar surface density profiles for the twelve reference models at the final time-step of each simulation. Plotted is the 
projected surface density a [g/cm 2 ] against the logarithmic radius [kpc]. Blue dots: stellar particles radial averages; dashed line: best 
Sersic fit (see text for details). The black diamonds are placed at the effective radii r e . The red squares are placed at the smallest value 
reached by adaptive softening lengths, thus indicating the effective resolution of the simulations. Left to right: high density, intermediate 
density, low density, very low density. Top row: high mass; intermediate row: intermediate mass; bottom row: low mass. 
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Figure 24. Density profiles for the twelve reference models at z ~ 1. Plotted is the density p [g/cm 3 ] against the logarithmic radius [kpc]; 
blue dots: stellar particles; cyan dots: DM particles; dashed lines: best Hernquist and NFW fits. Left to right: high density, intermediate 
density, low density, very low density. Top row: high mass; intermediate row: intermediate mass; bottom row: low mass. 



which gas collapses and forms stellar systems. Historically, 
there are two possible and competing theoretical scenarios 
to describe the assembly of the stellar mass of an elliptical 
galaxy (excluding the possibility of merger events). In the 
first one, stars originally form far away from the effective 
radius of the galaxy and subsequently accrete onto its inner 
regions, in a non-dissipative fashion. In the other case, first 
the gas flows into the central regions and then it is turned 



into stars. If the correct scenario is the first one, smoothing 
of the central density cusp of the profiles is expected, because 
the infalling stellar clumps loose orbital energy through dy - 
namical friction and the halo heats up (|Bertin et al.ll200a ). 
In the second case, the central distribution is expected to 
steepen up, because the gas radiates away it s internal en- 
ergy and an adiabatic contraction takes place l|Gnedin et al.l 
2004). Obviously, both processes are expected to play a role 
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in reality, and this is indeed what we see in our models, 
in which however the final density profiles show a plateau 
in the central regions, as expected in the "infalling stars" 
model. 

Our results also clearly suggest that if stars are formed 
too early, the resulting stellar systems are too shallow, end- 
ing up in "wrong" profiles: stars will tend to conserve their 
velocity dispersion as they fall into the potential wells and 
will therefore have elongated orbits, which is not the case 
when gas is let dissipate part of its kinetic energy and reach 
very high densities before turning into stars. 

In this context, it is worth recalling that there is some 
observational evidence for a deficit in the central luminosity 
of a few early type galaxies as recently found in observa- 
tional data bv lCote et al.l (|2007l .their Fig. 1), and references 
therein). However, the observed flattening begins at smaller 
radii (a few percent of r e ), and it is present only in bright 
(i.e. massive) systems, whereas the small systems seem to 
present a luminosity excess. Indeed, looking at the panels in 
the bottom row of Fig. 1231 one can notice that the central 
profiles of low mass models show some scatter and tend to 
steepen in the very central regions (especially in the VLDLM 
model). Recalling that the observational information stands 
on the luminosity profile, given the presence of young and 
hence luminous stellar populations in the cores of our low 
mass galaxy models, (and their absence in the cores of the 
massive ones), we can speculate that the trends of the ob- 
served curves should be fairly reproduced by translating the 
mass profiles into luminosity into luminosity profiles. This 
is an interesting point, because the excess or deficit of light 
could be explained just in terms of the stars and DM orbits, 
without the presence of a super-massive black hole affecting 
the global dynamics in the center of the systems, as instead 
suggested in the studies quoted above. 
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Figure 25. Mass-radius relation, i.e. the relation log(r 1 / 2 ) versus 
log M12, where T\/2 is the half mass (light) radius, and logMi2 is 
the total stellar mass in units of 10 12 Mq. The contour lines dis- 
play the regions of the MR-plane populated by the same number 
of objects (the bar yields the correspondence between color-code 
and number of galaxies). The data are from the HB sample of 
iBernardi et alj l|2010f ). The green solid line is their linear best fit, 
log r i/2 = 0.54 log M s — 5 .25. The black diamonds are the data 
from lBurstein et al which also include the dwarf galaxies. 

Large squares: the twelve major models at the last simulated age; 
the same colors refer to same initial densities. Black: high density. 
Blue: intermediate density. Magenta: low density. Red: very low 
density. The thin dashed black lines are fits to the models (see 
text). 



5 THE MASS-RADIUS RELATION 

To conclude our analysis, we briefly consider the Mass- 
Radius Relation (MRR) for ETGs. An overview of recent re- 
sults concerning estimates of the star mass and dimensions of 
ETGs is given by the MRR sho wn in Fig. | 25 | which dis plays 
the sample of ETGs collected bv lBernardi et al.l l|20ld , green 
dots). The sample contains about 65,000 objects. We have 
also added the small group of galaxies (about 60 objects in 
total) studied long ago bv lBurstein et all l| 19971 ) to include 
also some dwarf galaxies. Owing to the huge difference in 
the total populations of the two sample cno cross checking 
of the data has been applied for the galaxies in common. 
This plot is only meant to provide an idea of position of real 
galaxies onto the MR-plane. We also plot the position of the 
reference models coding them with different colors according 
to the initial over-density. 

At a first glance, one would say that high and intermedi- 
ate mass models fairly agree with the observations, whereas 
the low mass ones apparently have too large radii with re- 
spect to their masses. Therefore if on one hand the dimen- 
sions of intermediate and high mass models agree with the 
data, thus improving upon the situation with the hierar- 
chical scheme that predicts too small radii, on the other 
hand because of deviation shown by the low mass models, 
the slope of the observational MRR cannot be accounted 



for. However, the t horough analysis of this issue made by 
ICh iosi et all |2012a) has shown that the observational MRR 
is the result of a subtle interaction among different concur- 
ring physical causes. 

We begin by noting that as expected models with equal 
initial over-densities fall along the straight lines logri/2 = 
a log M a + f3, with slope a ranging from 0.2 to 0.3 and /3 de- 
pending on the redshift ft oc (1 + z)~ l . Slopes close to 0.33 
are indeed what is predict ed for collapsin g spherical per- 
turbations of DM (see e.g. iFan et aLl l2010). The deviation 
from the pure logri/2 = alogM s + P relationship is larger 
towards galaxies of lower mass an d it is due to the c o mplex 
bar yon physics; the disc ussion by IChiosi fe Carrarol (|2002f ) 



and 



Chiosi et al 



Chiosi et al 



(2012a) on the issue still applies also here. 
i|2012al ) argue that the observed MRR for 
ETGs is the result of two complementary mechanisms. On 
one hand, the mass function of DM haloes hosting the vis- 
ible galaxies gives (i) the typical cut-off mass at which, at 
any redshift, haloes become "common" on a chosen spatial 
scale, and (ii) the typical epoch at which low mass haloes 
begin to vanish at a rate higher t han that at which they 
are born because of merger events iLukic et all (|2007l ). On 
the other hand, these constraints define two loci (curves) on 
the MR-plane, because to each mass and redshift a typical 
dimension (i.e., radius) can be associated (using the basic 
relation between mass and radius of a collapsing object). 
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If the typical dimension of a galaxy is somehow related to 
that of the hosting DM halo (as our NB-TSPH models seem 
to suggest), then the region of the MR-plane between the 
two limits fixed by the halo mass function is populated by 
galaxies whose dimensions are fixed at the epoch of forma- 
tion, and only those objects that are "possible" at any given 
epoch may exist, populating a narrow region of the MR- 
plane. Al l details of this rath er complicate analysis can be 
found in IChiosi et al.1 (|2012ah to whom the reader should 
refer. 



6 SUMMARY AND CONCLUSIONS 

We have presented a numerical investigation on the star for- 
mation histories of ETGs, performed with the aid of hy- 
drodynamical NB-TSPH simulations of galaxy formation in 
cosmological context. By varying the initial properties of a 
single model in terms of total mass and initial over-density, 
we investigated how the initial configuration influences the 
formation of the system along the Hubble time. 

The main result of our study is that the star for- 
mation history of a galaxy is directly dependent on its 
initial total mass and initial over-density. Massive galax- 
ies (Mtot — 10 13 Mq) build their stellar content passing 
through an initial burst of activity at early times. In con- 
trast, low mass systems (Mtot — 10 9 Mq) continue to form 
stars throughout their lifetime, with episodes of intense ac- 
tivity separated by period of relative quiescence. Moreover, 
galaxies with the same initial mass have different star for- 
mation histories depending on their initial over-density: the 
deeper the perturbation, the more peaked, early and intense 
is the star forming activity. This is particularly true for inter- 
mediate mass galaxies (M to t — lO 11 Af0). Thus, denser ob- 
jects, which may be identified with cluster and group central 
galaxies, have more peaked and earlier burst of star forma- 
tion, while less dense objects (i.e., field and satellite galaxies) 
usually have delayed and lower ac tivity. This is larg e ly con - 
sistent with previous anal yses by Chiosi fc Carrarol (2002) , 
IClemens et~atl (120061 ). and lThomas et al.l (|2005h to mention 



a few. 

The structural properties of our models are in good 
agreement with those of the observed real galaxies, with 
some important differences (e.g. the density profiles in the 
central regions, the mass-to-metallicity relation) which we 
ascribe to the unusual choice we have made for the star for- 
mation dimensionless efficiency parameter, esF = 1. 

The picture emerging from this study is that an early- 
hierarchical, monolithic-like intense activity of star forma- 
tion at high redshift is possible starting from physically 
sounded cosmological conditions, and is fully adequate to 
explain the observed features of massive elliptical galaxies. 

Major mergers are certainly possible, but they are not 
strictly required. It is indeed easy to understand that on the 
contrary they would alter the subtle game among different 
physical ingredients all contributing to shape the body of 
properties and their their mutual tight correlations shown 
by early type galaxies such as for instance th e MRR it- 
self (IChiosi fc Carraroll2002l;IChiosi et alj|2012bl ). Along this 
line of thought iBuonomo et al.l l|200ol ) showed how a major 
merger would shift the position of the parent objects to- 
wards the upper region of the MR plane, where no observed 



galaxy is found. Moreover, chemical features of the elliptical 
galaxies cannot be reproduced if (wet) merging events are 
considered. F i nally, the shape of the galaxy ma ss function 
|Bundv et al.l 120061 ; IPerez-Gonzalez et al.1 120081 ) rules out 
the possibility that late major dry mergers are responsible 
for the mass assembly of massive systems, out of smaller but 
already passively evolving objects. Our models show instead 
that massive galaxies can form out of single, peaked pertur- 
bations, consistently with the cosmological background. 

Another interesting point to note is that the stellar feed- 
back seems to be sufficient to quench the star formation pro- 
cess in massive objects without invoking more exotic sources 
of energy, like Active Nuclei, provided it is modeled with suf- 
ficient accuracy. Since the models are evolved in void, the 
lack of late-time infall of gas could explain why the presence 
of an AGN is not strictly necessary. In relation to this, one 
could argue that strong galactic winds moving outward at 
high velocity are probably sufficient to inhibit flows of ma- 
terial towards the inner regions. However, the occurrence of 
Active Nuclei cannot be excluded. 

Furthermore, we are also able to recover the com- 
plex, episodic star formation histories typical of many dwarf 
galaxies. The galactic breathing phenomenon is a natural 
consequence of an accurate treatment of the stellar feedback 
process. 

We conclude that the fate of a proto-galaxy is essen- 
tially determined by its initial conditions in terms of mass 
and over-density, while external factors such as encounters, 
mergers and disruptions, while substantially altering the 
evolution of the involved systems, are not a fundamental 
ingredient in the global evolution of the early type galaxies 
population. 
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